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ABSTRACT 


A weather radar signal simulator that produees an output eonsisting of a veetor of 
I and Q values representing the radar return permits investigation of the performanee of 
different estimators for the weather signal parameters and their sensitivity when varying 
radar parameters and precipitation models. Although several empirical statistical models 
are available to describe precipitation behavior, the creation of a physical model enables 
adaptation to actual data (e.g. rain rate, wind shears) thereby making it possible to apply 
and examine different scanning schemes, especially rapid scanning schemes. A physical 
model allows gradual improvements to realism to study the effects on the radar return for 
different phenomena. A Weather Radar Signal Simulator has been developed in 
MATLAB. Several different functionalities have been implemented allowing for stepped 
frequency, multiple PRFs, pulse compression using a chirp, and variation of both weather 
and radar input parameters. Post processing capabilities include autocorrelation and FFT 
(for single PRF only); estimation of weather parameters such as reflectivity factor, Z; 
average doppler, radial velocity, and velocity spread; pedagogical plots including a 
Phaser plot of phase change over time and a velocity histogram, instantaneous observed 
reflectivity and power for each pulse over time. 
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EXECUTIVE SUMMARY 


Electromagnetic emission is increasing in all surroundings not only in the military 
arena making it necessary to understand, interpret, and sometimes filter information 
created by man or nature. In an attempt to learn more in one research area, studying 
other areas might provide the insight needed to break new ground. Electromagnetic 
returns from weather are in one area of radar science regarded as clutter while in another 
it as a way of gaining information and making estimates and predictions. Exploring the 
weather side of radar science provides an insight into weather returns as source of 
information and also improves knowledge about effects that can lead to a better 
understanding of how to overcome unwanted weather returns. 

Creating a weather radar signal simulator that produces an output consisting of a 
vector of I and Q values representing the radar return, permits investigation of the 
performance of different estimators for the weather signal parameters and their sensitivity 
when varying radar parameters and precipitation models. In weather radar applications 
the return signal samples can be used to estimate reflectivity, velocity and velocity 
spread. These are the zeroth, first and second moments of the Doppler spectrum. 
Reflectivity can be used to estimate the rain rate of the specific volume return. A number 
of mathematical models relating rain rate to reflectivity exist and can be used when an 
estimated power return is available. There are also a number of measures and statistics 
for precipitation that have resulted in functions describing rain drop size distribution for a 
certain rain rate, terminal velocities of rain drops, etc. 

Eimiting the amount of samples required to accomplish small variance when 
estimating average power calls for independent samples meaning that the return samples 
are de-correlated. On the other hand, measuring Doppler calls for coherent, correlated 
samples, this typically requires a higher pulse repetition frequency (PRE). The velocity 
spread effects the weather signal correlation time and therefore drives sample time 
required to arrive at an estimate with specific variance. Although several empirical 
statistical models are available to describe precipitation behavior, the creation of a 
physical model enables adaptation to actual data (e.g. rain rate, wind shears) thereby 

xvii 



making it possible to apply and examine different seanning sehemes, espeeially rapid 
seanning sehemes, as well as to ereate output graphieal representations that ean serve a 
pedagogieal purpose. 

This research was undertaken to build a weather radar signal simulator delivering 
I and Q values from a physical representation of rain drops, which has been 
accomplished. A previously developed simulator was initially evaluated and verified 
serving as the basic building block when developing the final version. The Radiation 
Integrals were used to derive the electromagnetic scattering from rain drops 
approximated as dielectric spheres. Without accounting for multiple scattering, 
attenuation, or coupling effects, the scattered electric field from each drop is summed for 
every pulse providing the RCS of the rain in the radar resolution cell as well as a complex 
voltage return representation. Using the Marshall-Palmer drop size distribution, an initial 
position for every drop is determined within an extended radar resolution volume 
enabling the drops to fall through the radar beam for all samples. Between pulses the 
drops are moved according to their respective terminal velocities based on the Atlas- 
Ulbrich approximation and wind. To allow for width in the velocity spectrum, a 
Gaussian spread is applied to drop velocities that determine how the average velocity is 
spread. Also a Gaussian spread, that is size dependent, is applied to the velocity and it 
represents the spread of velocity due to difference in size. The simulator power return is 
weighted using the Radar Range Equation for point targets, considering all drops as point 
targets within the resolution cell. The zeroth moment estimate is derived by averaging N 
power samples. 

Several different functionalities have been implemented allowing for stepped 
frequencies, multiple PRFs, pulse compression using a chirp, and varying input 
parameters. Post processing capabilities include autocorrelation and EFT (only for single 
PRF); evaluation of weather parameter estimators such as average reflectivity factor, Z; 
average doppler, radial velocity, and velocity spread. Pedagogical plots including a 
Phaser plot of phase change over time and a velocity histogram, instantaneous observed 
reflectivity and power for each pulse over time, have also been implemented. 



I. 


INTRODUCTION 


A. BACKGROUND 

Electromagnetic emission is increasing in all surroundings not only in the military 
arena making it necessary to understand, interpret, and sometimes filter information 
created by man or nature. In an attempt to learn more in one research area, studying 
other areas might provide the insight needed to break new ground. Electromagnetic 
returns from weather are in one area of radar science regarded as clutter while others use 
it as a way of gaining information and make estimates and predictions. Exploring the 
weather side of radar science will not only provide an insight into weather returns as 
source of information but also improve knowledge about effects that can lead to a better 
understanding of how to overcome unwanted weather returns. 

Creating a weather radar signal simulator that produces an output consisting of a 
vector of I and Q values representing the radar return, will permit investigation of the 
performance of different estimators for the weather signal parameters and their sensitivity 
when varying radar parameters and precipitation models. In weather radar applications 
the return signal samples can be used to estimate reflectivity, velocity and velocity 
spread. These are the zeroth, first and second moments of the Doppler spectrum. 
Reflectivity can be used to estimate the rain rate of the specific volume return. A number 
of mathematical models relating rain rate to reflectivity exist and can be used when an 
estimated power return is available. There are also a number of measures and statistics 
for precipitation that have resulted in functions describing rain drop size distribution for a 
certain rain rate, terminal velocities of rain drops, etc. 

To minimize the number of samples required to realize small variance when 
estimating average power independent samples are required meaning that the return 
samples are de-correlated. On the other hand, measuring Doppler calls for coherent, 
correlated samples, this typically requires a higher pulse repetition frequency (PRF). The 
velocity spread effects the weather signal correlation time and therefore drives sample 
rate. Although several empirical statistical models are available to describe precipitation 
behavior, the creation of a physical model enables adaptation to actual data (e.g. rain rate. 
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wind shears) thereby making it possible to apply and examine different seanning 
sehemes, espeeially rapid seanning sehemes, as well as to ereate output graphieal 
representations that can serve a pedagogical purpose. The output signal should therefore 
correspond to the values expected using the empirical models. The simulator should also 
enable input of data to create variation in reflectivity, velocity spread and Doppler. 

The simulator is based on the reflection of each drop creating a scattered electric 
field. A Marshall-Palmer drop size distribution is used. The sum of the fields scattered 
by all drops at any particular instant of time gives the instantaneous scattered field, total 
radar cross-section and instantaneous received power. Motion is imparted to each drop 
using a terminal velocity based on drop size and wind field. The motion leads to a 
rearrangement of drops in the sample volume and thus to a new value of instantaneous 
received weather signal power at the next sample time. The effect of electromagnetic 
coupling among drops is investigated. 

B. OBJECTIVE 

The main objective was to build a weather radar signal simulator that produces an 
output consisting of a vector of I and Q values of the radar return. From this, average 
power, Doppler and Doppler spread can be derived. To guide the research and provide 
tools for verification and validation of the simulator the following questions were 
addressed; 

1. Is it possible to build a weather radar signal simulator based on a physical 
model of a spatial region containing raindrops? 

2. Can realistic motion that depends upon the type of weather system be 
imparted to the raindrops? 

3. Can the output of the weather radar simulator be used to study spatial and 
temporal sampling schemes for constant or stepped frequency sampling 
pulses? 
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4. Can the weather radar simulator be used to study the use of pulse 
compression and range averaging as a means of rapidly obtaining 
independent samples? 

5. Can the output of the weather radar signal simulator be used to produce 
displays of pedagogical interest? 

6. Can the output of the weather radar signal simulator be used to study the 
performance of estimators for the first three Doppler moments? 

7. Can the output of the weather radar signal simulator be used to study the 
utility of estimators for higher order Doppler moments? 

C. APPROACH 

The simulator was developed by building modules that represented a functionality 
or physical behavior, testing and verifying the module, and then placing the module into 
the final version of the simulator. Having a prototype from France provided valuable 
ideas and experience for developing the simulator. Initially this version was translated, 
evaluated, and respective function or module output was verified against empirical 
models based on statistics of real life measurements. When all improvements were 
implemented to the French version, new modules were created aiming at addressing the 
research questions posed. The final product was optimized to give accurate results 
structured and presented in a pedagogical manner. During the development Graphical 
User Interfaces (GUIs) were of secondary concern. 

When the simulator was working, i.e. it has passed the necessary testing for 
verification; the application phase began. Research questions three to seven were then 
addressed. Results were again compared with empirical models built on statistics of real 
life measurements. Questions not addressed were formulated for future investigation by 
others. 

Although the subject title implies that meteorology will be extensively covered, 
the goal of this work was to provide a tool to evaluate the effects on the radar output of 
different weather systems. Precipitation physics was used to build a physical model, but 
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initially weather system input was ineluded only as a eonceptual mathematieal veetor 
model to provide flexibility to later implement real models of weather. 

D. RELATED WORK 

M. Gosset, J. Nieol, and A. Sanchez submitted an abstract for The 6th Conference 
on Hydrological Applications of Weather Radar describing a simulator with similar focus 
[1]. After personal contact with M. Gosset we learned that the simulator was unfinished 
and later the MATLAB code was sent to us to develop at NFS in cooperation with M. 
Gosset. 

J.S. Marshall and Walter Hitschfeld [2] thoroughly investigated theoretical effects 
such as fluctuation echoes from randomly distributed scatterers providing mathematical 
proofs and suggestions for pulsing and scanning schemes to improve estimate accuracy. 
The paper defines probability distributions useful to interpret fluctuating returns and 
describes different methods to achieve independent data using difference in time, range, 
and frequency, to mention a few. 

Richard J. Doviak and Dusan S. Zmic [3] provide extensive guidance in the 
weather radar area presenting tools and explanations of all aspects investigated in this 
thesis. Most background data has originated from thoughts, theories, and descriptions 
provided by their fantastic book. Building a simulator calls for in depth interpretation of 
radar hardware solutions, which in large part can be found in their work. 

Finally, John B Sandifer [4] developed scanning strategies for research and 
operational applications. Interesting temporal schemes were presented aiming to provide 
methods to rapidly scan and update large volumes using Phased Array Doppler radars. 
The use of stepped frequency to acquire independent samples without having to wait for 
the weather to reshuffle enough for de-correlation is explored and presented as one 
successful way of achieving a fast scanning capability. 
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E. THESIS OUTLINE 


The thesis is organized as follows: 

Chapter II, titled Weather Radar and Weather/Rain theory, introduces the basic 
principles of radar and the specific functions in weather applications. The physics of 
precipitation is summarized with emphasis on electromagnetic effects (scatter), motion 
behavior, and drop size distribution. The chapter also describes Radar Cross Section 
(RCS) of rain. The basis of the simulator is the radar signal return from the rain, which is 
why the electromagnetic effects will be examined to provide the foundation for the 
simulator design. 

Chapter III covers simulator design considerations and the building of the 
simulator. In this chapter the scope and the limitations of the simulator are presented 
based upon the research questions that need to be addressed and the abilities of the 
simulator building tool of choice (MATLAB). The chapter further describes the building 
of the simulator. This chapter covers both the development and testing of the simulator. 
MATLAB is the computer language of choice and the simulator is initially developed in 
modules, to simplify development and testing. Ending the chapter the simulator parts are 
described. 

Chapter IV provides analysis of applications. The final test or validation of the 
simulator where the research questions 3-7 are addressed. 

Chapter V presents the results from both the development and the applications 
covered in the thesis. 


Chapter VI contains conclusions and recommendations for future work 

Appendix A Derivation of RCS for single drop 

Appendix B Microwave Studio simulator setup 

Appendix C Simulator setup for the Weather Radar Simulator 

Appendix D MATLAB code 


5 



THIS PAGE INTENTIONALLY LEET BLANK 


6 



11. WEATHER RADAR AND WEATHER/RAIN THEORY 


Laying the foundation for the simulator, the basie theory is eovered in this ehapter 
ending with a statement of ground truth for later eomparison and verifieation. 

A. RADAR PRINCIPLES 

Radar or RAdio Deteeting And Ranging uses the return of transmitted 
eleetromagnetie energy from a transmitter to deteet and loeate objeets of interest. 
Although the basie prineiples of radar have been known for deeades new applieations are 
developed eonstantly as a eonsequenee of the exploration of different areas of interest. 
The basic function of a radar is provided by a transmitter that generates an 
electromagnetic signal into space by the means of an antenna. Objects in the signal path 
scatter or reflect part of the power projected by the transmitter back towards the receiver. 
By measuring the return in terms of angle of incidence and time to the return, power 
levels, etc. different factors can be estimated depending on the application. [5] 

To be able to detect movement or to calculate average power, more than one 
return pulse must be measured leading to the idea of Pulse Repetition Frequency (PRF). 
By selecting different PRFs some quantities will be become ambiguous while others 
become the unambiguous. 

I. Estimation of Range 

Finding the range to the object of interest is merely a question of measuring the 
time it takes for the signal to travel from the transmitter to the target and back to the 
receiver, or more precisely from and to the antenna. Knowing that the electromagnetic 
energy travels at the speed of light, c, which in free space is approximately 3x10* m/s, 
the distance R is given by 

R = f ( 2 . 1 ) 

where T is the time from the antenna to the object and back to the antenna. 
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To ensure unambiguous measurements of range the time differenee between two 
eonseeutive pulses must be at least the time it takes for the signal to travel to the range of 
interest and baek. The maximum ambiguous range, is thereby 


^ _cT^ep _ C 

“ 2 24 

where is the pulse repetition period and 4 the PRF. 


( 2 . 2 ) 


2 Estimation of Power Density 

a. Antenna Gain and Beam Shapes 

An important part of a radar system is the antenna. The antenna is the 
impedanee mateher making it possible to transmit and direct energy into space with 
minimal losses. Some antennas provide the ability to direct the energy, making the 
power density higher in a limited solid angle while being lower elsewhere. This ability, 
called directivity, is converted to a measure of power per unit solid angle radiated in a 
particular direction and is compared to the same measure but without directivity. The 
ratio is called the directive gain of an antenna and is a function of direction and it does 
not include losses in the antenna. 

Apart from resulting in a higher Gain, the ability to shape the beam can be 
helpful to enhance angular resolution, since it will limit the volume of illumination. An 
isotropic antenna has a power density at a range R of 


P = —^ 


AttR 


(2.3) 


where is the transmitted power, while an antenna with gain has the power density of 




47rR 


2 * 


(2.4) 


b. Radar Range Equation 

To be able to estimate the return power from an object, a parameter called 
Radar Cross Section (RCS) must be defined. The RCS can be defined as [6] 
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_ Power reflected to receiver per unit solid angle 
Incident power density/47i 

The scattered density power at the receiving antenna is given by 

P.G ( 


P 

ant 


[W/m^]. 


Accounting for the capture area of the radar antenna, , results in received power 

By using the relationship between transmit gain and receive effective area, 

4;rA„, 


(2.5) 


( 2 . 6 ) 


G = 




/I' 


(2.7) 


where X is the radar wavelength, received power can be expressed as 

P,G,GM" 

Gee , \3 4 ( 2 - 0 ) 

(4;r) R^ 

ignoring all system and propagation losses. Since losses are expected let us for now add a 
one-way propagation loss term L rendering 


p _ P,G,GM^ 


{Any R^L^ 


(2.9) 


3. Estimation of Velocity 

Pulsed radar can, by measuring the phase change between consecutive pulses, 
estimate radial velocity. The effect, called Doppler, is the frequency change that occurs 
due to the relative movement of a target with respect to the radar. Considering a two-way 
propagation path renders a total phase change of 

2 IP 

(l) = 2n -. (2.10) 

X 

In terms of angular frequency the result is 
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( 2 . 11 ) 


dt A dt X 

where is the radial velocity and is the doppler frequency [4], 


From a pulse pair, using a coherent detector, the phase change can be extracted 


using 


V = I + jO = A 

n n d n 

arg{y„} = ^„ 


V V = A A e 

''n+V n ^n+l^n'^ 




leading to the single pulse pair velocity estimate 


( 2 . 12 ) 


^ Ft/ t/*1 '^■^p Ft/ t/*1 

^ =-arg j = ^ arg |y„,jy„ j 


4>^r„ 


(2.13) 


where T is time between pulses and / is pulse repetition frequency. 


B. WEATHER RADAR 

The radar application that will be explored in this research is the ability to use the 
radar to estimate weather parameters such as rain reflectivity and wind speed. Most radar 
applications consider rain as clutter and find ways to filter out those effects. In weather 
applications the scatter from precipitation is used to compute the estimates mentioned 
earlier. 


1. Sample Correlation 

The return signal power level from weather will vary over time. To 
compute an estimate of reflectivity with small variance a number of independent samples 
are summed and averaged which means that the target volume must reshuffle enough to 
ensure de-correlation during the observation time. The power spectral density of a 
meteorological signal is approximately Gaussian [7] and can be written as 
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lal 

5(/) = V^ J (2.14) 

where / is the doppler frequeney and cr^ is the doppler speetrum varianee. Taking the 
Fourier transform of (2.14) and normalizing results in the eorrelation eoefficient 




r 

2 ( 7 ^ 


(2.15) 


where is the time speetrum varianee and is related to aj by 


and 


1 

27t(Jj 



(2.16) 


(2.17) 


By using the eorrelation eoeffieient a measure of minimum differenee in time ean 
be estimated ensuring dependent or independent samples. Nathanson elaims that for 
independent sampling /7 (t)<0.02 is required whereas /7 (t)> 0.15 will insure 
dependenee [7]. Doviak and Zmic [3] set a eorrelation threshold for eoherenee at 
(r) > e ' whieh eorresponds to a value of p{t) > 0.6 . Using this threshold renders 


r^<0.25— (2.18) 

KO^ 


for highly eorrelated samples. For independent samples using [7] > 0.70- is 

ncj^ 

required. 


2. Depth of Volume 

Estimating range to a volume target, eompared to a point target, still involves 
measuring travel time. The range to every raindrop will not be determined sinee the 
return signal will be the sum of several seatterers. Measuring a volume takes pulsewidth 
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into account as one of the parameters determining its size. A pulsewidth of T renders a 
depth of the return volume of 


dy 


CT 


(2.19) 


3. Estimation of Power Density and Reflectivity Z 

a. Antenna Gain and Beam Effective Solid Angle 

The direetivity of an antenna is a funetion of both azimuth and elevation 
angles. The radar beam has a width that has to be speeified to enable an estimation of the 
magnitude of the volume return. Estimates of a point target involve the angle off bore 
sight to get the aeeurate gain level, while a volume target, like rain, takes into aeeount the 
total gain over the volume of interest. Estimation of the gain for a volume target uses the 
antenna pattern funetion /(^,^) integrated over the solid angle. Sinee the radar 

equation is a power relation the pattern funetion needs to be squared, also, due to the two 
way propagation, the power pattern function is squared again leading to 

( 2 . 20 ) 


Using a Gaussian Pattern Model provides a good approximation of the 
main beam between the 3 dB points. Extending it to both azimuth and elevation renders 

( 2 . 21 ) 

where ^— and = ——, and 0 ^, are the half power beamwiths in the 

4 In 2 4 In 2 

orthogonal planes. Performing the integration over the solid angle renders [8] 





If' 


] 


dOd(l)- 


Att 


81n2 


(1-^) 


( 2 . 22 ) 


where ^ < 0.034 and represents the eontribution from the sidelobes. It eorresponds to 
0.15 dB error allowing for an approximation ignoring the term (l-^) rendering 
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(2.23) 




7t0^(p^ 

81n2’ 


h. Weather Radar Range Equation and Estimation of Reflectivity 
To define the weather radar range equation we return to the original 

formula 




(2.24) 


In the weather radar equation RCS (<t) is replaeed by RCS density rj, whieh is the 
expeeted radar eross seetion density per unit volume. Sinee the target is a volume target, 
the expeeted return needs to be integrated over the whole volume, whieh ean be divided 
into the solid angle and depth of the return. Letting the power level at the target be 


target 


P,G, 

47rR^ 


dA- 


P.G. 




‘-dQ. 


and integrating the return power over the return volume 


renders 


R+— 

3 2yr-, ^ 2 


n -iiP {-dR P,X‘G,G,cWA^ p.,,. 


R n 


CT 


reeognizing that R » —. 


Radar meteorologists use refleetivity faetor, Z, instead of instead of RCS 


density where 


n=YKtz. 


( 2 . 26 ) 


and 


AV^ 


(2.27) 


i&V 


Substituting (2.26) in (2.25) and solving for Z renders the reflectivity estimator 

1024RMn2/l" 


1024;r'RMn2 Z' 


rec „ Q 2 


= P ■ 

1 2 rec 


P,X^Gfiacre,(p, fG,G^cT0,</)y 


(2.28) 
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For an assumed continuous distribution of drops sizes with diameter density 
N[D) (2.27) becomes 

Z = N {D) D^dD. (2.29) 

Solving the integral form for an assumed Marshall-Palmer Drop Size Distribution, 
DSD (Equation (2.35)), yields 

Z = ]D‘N{D)dD= (2.30) 

0 (4.1i?-°-'') 


Drops do not exist in sizes ranging up to infinite diameter (see later sections for 
details) but are likely to be limited by a maximum diameter. Doviak and Zrnic [3] state 
that Equation (2.30) overestimates Z and suggest a truncated DSD with D = . In this 

case 


Z = 




(4.1i? 


-Z(7,a) 


(2.31) 


where a = KD^^ and Y{l,a) is the incomplete Gamma function. Equations (2.30) and 

(2.31) provide an analytical means for computing Z and providing a reference as ground 
truth reflectivity against which simulator estimates can be compared. 


4. Estimation of Velocity 

Measuring phase-change over time for a volume target demands coherence 
between pulse pairs. Phase change will also occur due to other effects rather than radial 
movement. To obtain an estimate of radial velocity with small variance, a number of 
sample pairs need to be summed and averaged. Building on (2.13) we obtain the multiple 
pulse pair estimator 


1 N-\ 

v-l^ 


n+1 


N-\ 


A^-1 


(2.32) 
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5. Estimation of Velocity Spread 

An estimate of the veloeity spread ean be computed using the velocity samples 
obtained from a pulse pair comparison (Equation (2.13)). The estimator for sample 
variance is given by 



and the estimator for standard deviation is 



(2.33) 


(2.34) 


C. PHYSICS OF PRECIPITATION 

Why and how raindrops form and what constitutes their shape and behavior is 
beyond the scope of this thesis, however, some basics in precipitation physics are needed 
to develop a physical representation of rain for the simulator. 

1. Size and Shape 

All drops examined are assumed to be falling at or close to their terminal velocity. 
Small drops, with a diameter D < 0.35 mm, are spherical in shape. Drops ranging from 
0.35 <D <4 mm have progressively flattened bases, and can be approximated by 
spheroids [2]. Larger raindrops tend to break up after collision and drops larger than 
about 10 mm in diameter are unstable and break up even without collision [9]. Also 
drops can cluster which will affect the drop size distribution. This phenomenon will not 
be covered in this thesis. For further reading [10] is suggested. 

2. Drop Size Distribution, DSD 

RCS density is the average radar cross section density per unit volume, as 
mentioned previously, and for independent scattering is the sum of the RCS of all 
scatterers in the return volume divided by the volume. The drop size distribution 
represents the diameter density (mm/m ) of drops for all diameters. The Marshall-Palmer 
drop size distribution [11], derived empirically, provides a general model, 
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(2.35) 


N{D) = NoQxp{-AD), 

A = 4.1i? 

Nq = SxlO^m^^mm 


Applying truncation as suggested in [3] with a maximum drop diameter of 

yields 


N{D) 


iN^cxpi-AD), 

lo^ 


D<D 

max 

D_<D' 


(2.36) 


3. Terminal Velocities 

The terminal veloeity of a water drop is dependent on the size of the drop. Atlas 
and Ulbrieh [12] derived an expression that relates terminal veloeity to the diameter of 
the water drop as 

v{D) = 1767D1'’’[cm/5] - 3S6.6D^f[m/s] (2.37) 

whieh is valid for 5x10 < D < 5x10"^ [m]. The result has been verified against 
measured data [3]. Fall speed is also dependent on air pressure but that is a seeond order 
effeet not eonsidered here. 

4. Velocity Spectrum Width 

There are several meehanisms eontributing to the spread of the veloeity speetrum 
both from nature and the radar itself. The radar antenna motion, turbulenee, differenees 
in fall speed due to drop size and wind shears all eontribute to the spread of the radial 
veloeity of preeipitation. Physieally modeling all eontributors demands great insight into 
both behavior and effects, and is beyond the seope of this thesis. However, to control de- 
correlation time a spread must be introdueed in the simulator. This is aeeomplished by 
imparting a random Gaussian distributed veloeity eomponent to drops. 

5. Rainfall Rate 

Rainfall rate determines depth of accumulated water per unit time, and ean be 
derived from water eontent and fall speeds. It ean be shown that the a eloud’s water 
density is 
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(2.38) 


M = '^p„\D^N{D)dD 

where is the water density of a drop. 

Turning this into rainfall rate renders 

R = -]D^N{D)v{D)dD. (2.39) 

6 0 


D. RADAR CROSS SECTION (RCS) OF RAIN 


I. Cross Section 

A eoherent measure of radar return ean be found by summing the eleetrie field 
eontribution from eaeh point seatterer in the return volume of interest. Raindrops ean be 
approximated as small dieleetrie spheres. It ean be shown (see appendix A) that the 
backseatter field of a small sphere, when the radius is mueh smaller than a wavelength 
(/I), is given by 


^ ^ li- P n6 


(2.40) 


where D is the drop diameter, 


(^r+2) 


and is the relative permittivity. 


This approximation is valid as long as the wavelength is mueh larger than the 
drop eireumferenee and the drop is of spherieal shape and is referred to as the Rayleigh 
approximation [3]. As deseribed under Physies of Preeipitation, drops tend to flatten 
when between 280 and 1000 pm whieh render a polarization dependent return. This ean 
be used to separate different sizes of precipitation when a dual-polarization radar is used 
to obtain the ratio of two orthogonally polarized returns [3]. This will not be included in 
the initial simulator work. 

To verify the return predicted by the simulator, the RCS model needs to be 
accurate. Not only does the model need to be accurate in terms of magnitude, the 
fluctuation over time must also be accurate to assure the model will produce the proper 
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variation matching the theoretical density function. For the magnitude test for accuraey, 
Microwave Studio has been used to predict the scatter from spheres with properties of 
raindrops. Microwave Studio uses the Finite Integration Technique to calculate RCS. In 
Microwave Studio a sphere was created with a diameter of 4 mm. Conductivity was set 
to 0.0001 S/m (pure water) [13]. The set up for all simulations was in most aspects the 
same and is described in Figure 1 and in detail in appendix B. 


z 



Figure 1. Simulation setup using Microwave Studio, single drop. k. is 
the incident field unit vector and E is the direction of the electric field. 

The calculations were made assuming an X-band radar. Although drops with 
D > 4 mm flatten out when falling, a spherical approximation was used allowing the 
Rayleigh approximation for the backscatter. For an X-band radar with frequency 

I |2 

9.4 GHz the theoretical backscatter of a single drop with \KJ^ = 0.93 , and a diameter of 
4 mm is 

(9.4x10’)' 0.93m^ 

= 1.124x10“'’ -^-59.5dBsm. 

The result using Microwave Studio with the same input values used in the above 
is shown in Figure 2 and Figure 3. 
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Farfield 'farlield (f=9.4) (pw]' RCS.AbsrTheta) 
90 



Figure 2. Farfield Scatter from a sphere of 4 mm in diameter in the 9 
plane for ^ = 90° and ^z) = 270°. 


Farfield 'farfield (N9.4) [pw]' RCS.AbslTheta) 



Figure 3. Farfield Scatter from a sphere of 4 mm in diameter in the 6 
plane for ^ = 0 and ^ = 180°. 
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In Figure 2 the level of the farfield scatter is equal in all directions, which is 
expected considering the electric field direction (x-direction) and the view of observation. 
The pattern is similar to a dipole where the farfield scatter has its maximum in the 
orthogonal direction from the electric field [15]. In Figure 3 the zero directions are inline 
with the polarization of the incident filed, which occurs when ^ = 0 and ^ = 180° with 
0 = 90°. 


2. Array Factor 

Multiple scatterers can be modeled as an array of scatterers even though they are 
not identical. One must realize that calculating the total coherent RCS of multiple 
scatterers is not a matter of summing the RCSs since RCS is acquired by summing the 
scattered electric field which then is squared to get the RCS. A simple approach is to 
calculate the electric field return from every scatterer and then add the contributions 
vectorally [6]. The phase change due to difference in location using a far field 
approximation is given by 

, (2.41) 

where represents the path difference for every drop with reference to , Zq„ and 

is given by (see Figure 4) 

Sa„ =Xo„sin6>cos^+yo„sin6>sin^+Zo„cos6>. (2.42) 


20 




X 

Figure 4. Far-field path difference due to shift from origin. 
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Figure 5. Simulation setup using Microwave Studio, two drops, k. is the 
incident field unit vector and E is the direction of the electric field. 

The value is ealeulated only for a direet refleotion not taking into aceount the 
cos^ 6 dependence due to angle of incidence. The two drop setup using D-A mm and 
d = A, rendered results presented in Figure 6 and Figure 7. In Figure 6 the farfield scatter 
pattern show a maximum in the orthogonal direction from the electric field, as expected 
since the electric field was along the y-axis (^ = 90° and ^ = 270° with ^ = 90°). In 
Figure 7 zero directions show for ^ = 0 and ^ = 180° with 0 = 30° and 0 = 150°, which is 
expected since 


sm[Nkd sin0) 

= 

sin 

\ 

2 - —/lsin0 
X j 


sin(4.;rsin0) 

sin [kd sin0) 

sin 

f 2;r , . y 

y X J 


sin(2.;rsin0) 


So for 0 = 30° and 0 = 150° the expression yields zero. Comparing the calculations, 
Equation (2.45), with the results from the Microwave Studio simulations reveals a 
difference of only 0.1 dB, which can be the effects of mutual coupling. The max values 
are called Bragg lobes and their locations are dependent on the ratio of d!X, their 
location with respect to each other, and the measuring point. 
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Fariield 'fariield (f=9.4| |pw]' RCS_Absn'heta) 



Figure 6. Farfield Scatter from two spheres of 4 mm in diameter in the 
-plane for ^ = 90° and ^ = 270°. 


Fariield 'fariield ((=g.4| |pw]' RCS_Absn'heta) 
90 



Figure 7. Farfield Scatter from two spheres of 4 mm in diameter in the 
0-plane for ^ = 0 and ^ = 180°. 
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The question of mutual eoupling was addressed using Mierowave Studio. The 
seattering pattern from the illumination of first two and then three spheres was examined, 
this time for different separation distanees (Figure 8). The values plotted are the results 
using the same drop setup as in Figure 5, i.e. D-A mm, 6 and ^ = 0 . The measures 
indieate eoupling effeets. However, the software is not optimized to ealculate low 
eonduetivity targets, as raindrops are. When the spaee between the drops inereases the 
program adds noise as the numerieal grid beeomes significant (see Appendix B). 
Although mutual coupling is likely to appear, available tools could not provide guidance 
how to provide realistic implementation. 


RCS for different distances between drops 



Figure 8. Monostatic Back Scatter of two and three drops of 4 mm 
diameter when spacing is varied from 0.1 A to 2A. 


3. Averaging over Time 

The weather radar return signal is composed of the scattered electric field from all 
drops in the resolution cell (see Figure 9). The total scattered field of the drops in the 
resolution cell is given by 
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(2.47) 


e,(o) = Se.(o) 

n-\ 

where N is the number of drops and O is the observation point. 


Rearranging, a = ^ renders ^ ^ , so summing the 

Ei Ei 

eleetrie field from all seatterers yields 

= (2.48) 


where y/ is given by y/^ = 2kgQ^ and is a funetion of r, 9 and (p (Equation (2.43)) [6] 
and a is a unit veetor giving the direetion of the seattered field. 


z 



X 


Figure 9. Radar resolution cell used in Simulator.! 

It ean be argued that y/ is a. random variable uniformly distributed on the 
interval [0,2.;r] thus, expeeted RCS will be 


1 Azimuth and Elevation angles are reversed compared to the mathematical conventions. 
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(2.49) 


N 

total ) “ 

n-l 

where ( ) denotes the expeeted value. So to get the average power estimate a number of 
independent samples need to be gathered and averaged. 

4. Absorption and Attenuation 

Apart from the fraetion of the ineident eleetrie field that rain seatters, rain also 
absorbs energy whieh eauses attenuation. This has greater effeet for higher frequeneies, 
whieh imposes limits on the possibilities of getting higher angular resolution by using 
short-wavelength radars. The effeets of attenuation ean be used in bi-statie setups 
estimating the weather parameters based on the level of attenuation. Bi-statie setups and 
effeets of attenuation will not be taken into aeeount in the initial simulator work. 

E. MATHEMATICAL MODELS 

This seetion deseribes the basie formulas that will be used to model the 
parameters required in the simulator to generate the return signal. The models for eaeh 
parameter are defined in Table 1. Although ehoosing Equation (2.30) as ground truth 
will probably ereate too high estimates of refleetivity [3], this will be the initial referenee 
as the simulator development starts. 
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Table 1. Mathematical models for Weather Estimators. 


Parameter 

Input or estimator 

Ground Truth 

Drop Size Distribution, 

DSD 

N{D) = NoQxp{-AD), 

A = 4.1i? ', 

= 8 X10^ m^^mm ' 


Reflectivity Factor, Z 

1024i?Mn2/l" 

“ P,G,G^cm^’\Kj 

Z = ]D‘N(D)dD= 

Average Radial Velocity 

— 1 p Ft/ t/*1 

Compare with input value 

Doppler Spectrum Width 

N , 

O-' = '=* , a = 4¥. 

'' A 

Compare with input value 
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III. SIMULATOR DESIGN CONSIDERATIONS AND 
DEVELOPMENT OF SIMULATOR 


Having described the required parameter models in the previous chapter the 
design and development of the simulator will commence. As stated, the simulator should 
allow for inputs of rain rate, standard radar values (e.g. frequency, azimuth and elevation 
angles, and pulsewidth), wind speeds and angles, and also velocity spread. Output signal 
parameter estimators of interest are average power, average Doppler velocity and 
Doppler spread. The output should preferably be in form of I and Q return voltage 
making it possible to post process the data to obtain the estimate of average power, 
Doppler velocity and Doppler spread. 

A. SCOPE 

To provide a tool to answer the research questions specified in the introduction 
chapter, the following is required: 

1. A physical model of a spatial region of raindrops must be created and later 
tested and verified against real life measures or existing models, 

2. A model for realistic motion of raindrops that are weather system 
dependent must be created where at least an interface between weather 
models and the simulator must be specified, 

3. A method of implementing spatial and temporal sampling schemes to 
study the output for constant or stepped frequency sampling, 

4. A method to study the use of pulse compression and range averaging as a 
means of rapidly obtaining independent samples, 

5. Displays of pedagogical interest to facilitate the presentation and testing of 
research issues mentioned, 

6 . A study of the performance of estimators for the first three Doppler 
moments, 

7. A study of the utility of estimators for higher order Doppler moments. 
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B. LIMITATIONS 


The initial design will not include; 

• Dual polarization 

• Absorption/Attenuation 

• Losses and noise 

• Mutual coupling and multiple reflections 


C. DEVELOPMENT 

Having a prototype from M. Cosset [1] provided valuable ideas and experience 
developing the simulator. The method used to develop the flnal version is described in 
Figure 10. 


Rewrite Modules 
& 

Construct Verification Programs 


French MATLAB code 
Module based (functions) 


4F0ffn)kf 


7JHfThlkhHKhJH 


Translation 
& 

Interpretation 


2 Gdira Qflhrmbd: 

7JHffh)khHiaiJH 


Simulation 

Speed 

& 

GUI 


A{D) = Vo exp(-AD), 
A = , 

Nq =8x 10’nr^ninr', 



Figure 10. Plan for development of Simulator using the French version. 
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1. Translation and Interpretation 

The French code consists of a main program that calls for functions that performs 
different tasks to support the main program. The functions are2: 

• Size 

• Approximation 

• Marshall-Palmer 

• Resolution Volume 

• Number of drops 

• Initial Position 

• Summation 

• Movement 

• Atlas Ulbrich 

The input parameters are shown in Table 2. 


Table 2. Simulator input parameters 


Radar Data 

Weather Parameters 

Transmit Power 

Drop size diameter limits 

Max Antenna Gain 

Drop Diameter resolution 

Range to resolution cell 

Rain rate 

One-way 3 dB Beamwidth 

Rain angles 

Pulsewidth 

Wind, speed and direction 

Frequency 


Pulse Repetition Frequency (PRF) 


Elevation and Azimuth 


Number of pulses for integration 



2 The code was translated by Professor Monique Fargues of NPS. 
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a. Function Description 

The main program starts by calling the Size-function that generates a 
discrete Diameter vector based on the diameter interval and resolution of interest. The 
function calls the Approximation-function that calculates the RCS for each diameter 
using Equation (2.40). The Marshall Palmer-function is called to calculate the number of 
drops for each diameter using Equation (2.35). The resolution is multiplied with the 
number of drops per diameter and stored in a matrix together with the Diameters and 
RCSs. The Size-function also has a part that uses a pre-set threshold to take away the 
smaller drops that account for 10 % of the total RCS. The matrix created in the Size- 
function consists of Diameter, RCS, Number of drops, total RCS per diameter, and 
Reflectivity Eactor structured as 


' A 



(T„ 

N{D,) 

... N{D„) 

a,-N{D,) 

■■ 

1 Z, 

J 


where the first values are set by the threshold and N by the upper Diameter limit. 


(3.1) 


The main program then calls the Resolution Volume-function to generate 
the resolution volume using the radar data. The function adds a margin to the resolution 
volume so that the drops can pass through the resolution volume as time passes between 
pulses. Both the borders for the increased resolution volume, henceforth called the Box, 
and the resolution volume are defined for later use. The main program uses the data to 
feed the Number of drops-function that calculates the number of drops considering the 
volume of the Box and changing the values in the Matrix. 


Using the Initial Position-function, the simulator uniformly distributes 
every drop inside the Box using the data from the Matrix. The positions are stored in 
files unique for each diameter using spherical coordinates to simplify comparison when 
determining whether a drop is inside the Resolution Volume or not. 


The main program then starts the simulation for the number of pulses 
specified. Eor all iterations the Summation-function calculates the return voltage and the 
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return power taking into account the position and RCS of each drop. The radar cross 
section is calculated summing the scattered electric field from each rain drop using 

■ (3.2) 

where is the sample time and s is the Diameter vector. 

The results are stored in different files depending on what data it is. 
Before next iteration the main program calls the Movement-Function to move the drops 
according to wind and fall velocities. To get the estimated terminal velocities for each 
diameter the Atlas Ulbrich-function [12] is used, which is based on Equation (2.37). 

The output of the simulator consists of four plots that show the Box with 
all initial drops plotted, the Resolution Volume with all positions the drops have during 
the simulation plotted, power return levels for every pulse, and the phase of every drop 
and every pulse. 

b. Verification 

To evaluate the accuracy of the simulator, verification programs were 
created to estimate errors of the modules of the simulator. The modules to examine were 

• Number of drops 

• Average power estimation 

• Doppler estimation 

• Doppler spread estimation 

(1) Number of drops. The original French program used a 
diameter interval between 0.05 mm and 5 mm with a resolution varying from 0.08 mm 
for the smaller drops (up to 1 mm) to 0.2 mm. Using the original rain rate set to 50 mm 
renders a histogram as seen in Figure 11. The threshold, that removes the smaller drops 
that account for 10 % of the total RCS, is also indicated in the same figure. 
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Drop Size [mm] 


Figure 11. Drop Size Distribution comparison for French simulator. The 
volume is set to 1 m^ and the rain rate to 50 mm. 

Eliminating the smallest drops will effeet the fluctuation of the 
return signal as well as the peak values, which in turn affects the average value. This is 
an effect of the array factor previously described that can make several small drops 
generate high return values when in phase. Also the choice of upper limits will affect the 
level of power return especially for higher rain rates, which will be examined in a later 
section. These values are plotted without making the number of drops to integer values. 
The simulator uses CEIL to round off the number of drops making the integer equal to or 
higher than the original value. Simulations show that the effect becomes significant for 
rain rates below 40 mm/hr, creating too high return power values (Eigure 12, RCS using 
original Erench simulator). To produce values closer to what the integral form generates 
using Equation (2.30), the span of diameters needs to be widened and the round off tool 
exchanged. 

(2) Average Power Estimation. To verify the average power 
return, a simplified simulator was built that didn’t include the movement of the drops. 
The movement was simulated by randomly re-plotting the drops for all iterations. Using 
the knowledge about the effects of choosing limits for the drops size, a simulation was 
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performed verifying the levels for average return power. Also the statistical distribution 
of the return signal was examined to assure it agreed with basic theory [1], [2]. The 
simulator was run for 1000 samples and for rain rates between 10 and 120 mm/h and a 
wavelength A-O.l m. 

The results (Figure 12) show that using the original French 
simulator makes the relation between RCS and rain rate almost linear while the new 
generator follows the integral form (Equation (2.30)). The new generator uses a wider 
span of diameters (up to 15 mm), does not apply a threshold and uses ROUND instead of 
CEIL making the value closer to the real value. 


RCS density of 1 of drops. 



Figure 12. Average RCS density comparison using different drop 
generators. 1000 samples were used and compared with the result using the integral 
form of reflectivity to calculate RCS density. A = 0.1 m. 


The new generator still produces values with errors up to 2 dB 
even with sample sizes as large as 1000 samples. To examine the cause of this error 
another version of a drop generator was constructed taking into account the error caused 
by the rounding off. The generator keeps track of the error and adds it at the end of the 

simulation. This also provided an idea of how to be able to improve the simulator speed 
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by reducing the number of drops while still retrieving accurate RCS and statistical 
behavior. Reducing the number of drops by dividing by a constant and after the 
simulations multiply it back into the calculations to retrieve reflectivity, number of drops, 
and rain rate is valid since all these operations are performed by integrating over the 
diameter. Treating the error like a constant results in a max error of about 0.4 dB, an 
error that seems to increase as the rain rate increases (see Figure 13). Since the simulator 
replaces the integration by a summation, errors will arise. Using Equation (2.30) renders 
a Z-value of 54.1 dBZ, for a rain rate of 100 mm/h. When using the same form of 
summation the simulator applies, the same estimation is 54.5 dBZ, which explains the 
error of 0.4 dB. 
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Figure 13. RCS density comparison using new drop generator with error 
correction and reduction factor. /l = 0.1 m. 

To get the average value of the return power, several independent 
samples are used. The RCS of a raindrop was derived using the radiation integral, which 
can be used when the return complex voltage is calculated. Since I and Q are random 
variables, independent and none of them dominant, the central limit theorem states that 
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their sum tends to be Gaussian, as long as the numbers of samples are large enough. The 
amplitude of the return voltage i.e. the absolute value of the I and Q voltage return, 

\y\ = [l" + Q"f" (3.3) 

is therefore expected to have a Rayleigh distribution [3]. Running the drop generator for 
10,000 samples renders results shown in Figure 14 and Figure 15. 


RCS density of InT^ of drops. Rainrate= 50. Number of drops= 4454. (^= 0, 0= 1.5708 

"56 I [ I I ^ I I ^ ^ 

- RCS density 

- Average RCS density 

- RCS density using eq (2.26) and Z=| D®N(D)dD 
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Figure 14. Plot of 10,000-sample-run of drop generator. 
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Figure 15. Rayleigh fit of data from 10,000-sample run of drop generator. 

In Figure 14 the average value of RCS is lower than the result 
using the integral form. This was explained earlier when the results from summing over 
all diameters was eompared with the integration from zero to infinity. The Rayleigh fit in 
Figure 15 follows the expeeted distribution. The results were retrieved by plotting a 
histogram of the 10,000-sample-run. Both the histogram and the fit were normalized to 
allow for comparison. 

(3) Doppler estimation. The module that moves the drops takes 
into account the terminal velocities, determined by Equation (2.37) with an associated fall 
angle, and wind speeds in x, y, and z-dimension. The spread of terminal velocities 
together with the fall angle will contribute to the spread of the Doppler spectrum. To 
evaluate the mean Doppler frequency a pulse pair algorithm is used to record the phase 
change over time using Equation (2.32) and comparing the phase extracted from two 
consecutive returns pulses. Applying this on the drop generator and the movement 
module enables an estimate of the mean frequency. A simulation was run for 30 pulses 
(29 pulse pairs) with a PRE of 2000 Hz. Wind speed was set to zero and the rain fall 
angle to 30° with respect to the z-axis (for detailed specifications see appendix C). The 
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simulation showed an average Doppler frequeney =-57.3 Hz or =-2.85 m/s in x- 

direetion (radial direetion) and standard deviation of (7^ = 26.8 Hz, whieh represents 

=1.34 m/s (see Figure 16). The spread was indueed by the differenees in terminal 

veloeities for different drops sizes. To estimate the aeeuraey, ealeulations were made 
using a 2 mm drop for whieh 

V (D) - 3S6.6D°f = 386.6 ■ 0.002°'" = 7.89[m/ 5 ]. 

Taking the fall angle into aeeount renders 

^x-dir =V'Sin(30°)eos(0) = -3.0l[m/5']. 

To further verify that the simulator generates aeeurate results, the 
fall angle was removed and only wind in the x-direetion was used. The same inputs were 

used adding wind of -5 m/s rendering fy =-105.1 Hz, whieh eorresponds to 
fy =-5.3 m/s as seen in Figure 17. The Doppler spread was (7^ =7.9 Hz (&^ = 0.4 

m/s). Note that the first speetrum (Figure 16) is skewed towards zero indieating that the 
smaller drops, eonstituting the majority of the drops, bias the result. Also the limited 
number of samples (29 pulse pairs) ean have effeets on the shape of the histogram. 
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Doppler Frequency Histogram with Gaussian Fit 
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Figure 16. Doppler Frequency Histogram applying terminal velocities and 
30° fall angle. 30 pulses and 29 pulse pairs. 

Doppler Frequency Histogram with Gaussian Fit 
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Figure 17. Doppler Frequency Histogram applying -5 m/s wind in -x- 
direction. 30 pulses and 29 pulse pairs. 
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(4) Doppler spread. The simulation run above showed a standard 
deviation of 1.18 m/s, which is a measure of how much the average radial velocity of all 
scatterers of each pulse return spread over all samples. The data was fitted to a normal 
density function, which is expected for a power spectrum [3], to explore if the data could 
be considered Gaussian. The second run fits better with the Gaussian plot, which can be 
explained by the randomness of the process different from the first simulation where the 
spread is distorted by the differences in drop sizes and therefore terminal velocities. To 
further examine the power spectrum new algorithms were implemented as described in 
the upcoming section. 

2. Adding New Modules 

The verification programs have been used to validate the simulator but there are 
still improvements that can be made to achieve better accuracy and realism. First to 
recapitulate what has been accomplished up until this point; 

A drop generator that provides the number of drops given by the 
distribution of interest based on rain rate. 

An algorithm that calculates the complex scattered voltage return based 
on the position of each drop within the radar resolution volume. 

A module that introduces movement to the drops concerning fall speed 
with an associated angle, and wind speed. 

An algorithm that calculates the doppler frequency of the moving drops 
based on phase change between pulse pairs. 

An error correcting algorithm that also can be used to reduce the number 
of drops used in a simulation to make it faster. 


a. Range Weighting 

To provide a more realistic level of the return signal value, a range 
weighting function, based on the radar range equation, was built and implemented. As 
mentioned earlier, effects of absorption and other losses will not be considered, leaving 
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only the effects of difference in range. To account for range effects, Equations (2.21), 
(2.24), and (2.48) were used rendering 


^ weigthed 
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where a is a unit vector giving the direction of the scattered field. For the total 
Resolution Cell 
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where y/^ = 2kgQ^ and is given by Equation (2.42). 


b. Correlation 

Another estimate of the doppler from frequency the return electric field 
can be obtained by taking the discrete Fourier transform (here EFT) of the signal samples 
using the complex voltage return. The Doppler Power Spectrum is given by 

S(f) = \2(ftTjM 

M-\ (3.6) 

m-0 

where M is the number of samples. An estimate can be retrieved by taking the Fourier 
transform of the autocorrelation function 


1 "H'h* 

^(0 = ^ S V(m)V(m + l). (3.7) 

M ^=0 
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Computing the autocorrelation function and the FFT for the output in the 

example on pg. 39 produces the same average frequency estimate, / =-105.1 Flz as 
before (see Figure 18). 

Power Spectral Estimate 
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Figure 18. Doppler spectrum using FFT and the Fourier transform of the 
autocorrelation function. (They overlap each other). Peak gives average frequency 

estimate / =-105.1 Hz. 

Both autocorrelation and FFT preserve amplitude information, which has 
consequences interpreting average Doppler using these functions. Larger drops will 
affect the power spectrum to a greater extent than smaller since drop RCS is proportional 
to . Thus, any difference in the radial velocity of large and small drops will result in 
the power spectrum being skewed towards the Doppler associated with the larger drops. 
Using the pulse pair algorithm (Equation (2.13)) removes the amplitude information 
producing a more “true” average doppler frequency or radial velocity. 

c. Frequency Stepping 

To accomplish rapid volumetric scanning there are several methods 
proposed. One method to achieve de-correlation is by shifting in frequency. A fast 
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scanning scheme to speed up the aequisition of independent samples was proposed by 
Doviak and Zrnic [3]. They suggested a method where a pulse pair is sent with one 
frequeney assuring dependent doppler measures, then the frequeney is ehanged a 
minimum of l/r, for de-eorrelation to send the next pulse pair. Implementing this 
seheme into the simulator requires modules stepping between more than one PRF also 
taking into aeeount the differenees in wavelengths when eomputing the meteorologieal 
signal parameter estimates. In Figure 19, a three frequeney stepping seheme is 
introdueed that transmits a pulse pair at frequeney 1, and listens for the time set by the 
Doppler Pulse Repetition Time, PRTi, before ehanging to frequeney 2. Before the 
seeond frequency can be transmitted the frequency oscillator needs a settling time, here 
called delay. After all three frequeneies have been transmitted the radar returns to 
frequency 1. The waiting time will be dependent on the de-correlation time for frequency 
1, which in this case is PRT 2 rendering a true time delay before resending frequeney 1 of 

tdelay = ) 

where N equals the number of frequencies used. 

Pulse 1 Pulse 2 

Frequency 1 


Frequency 2 


Frequency 3 


Figure 19. Frequency scheme using stepped frequencies to obtain de- 

correlation. 
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d. Pulse Compression 

A benefit of Pulse Compression lies in the eapability of averaging over 
wider range where the gain in range resolution provides independent samples that ean be 
used to estimate the average power for the volume of the wider pulse. Implementing 
ehirped pulse eompression in the simulator ealls for the ability to step in frequeney as 
well as range. A seheme was developed enabling stepped frequeney and range (see 
Figure 20), where the Resolution Cell is divided into smaller range bins defined by the 
sub pulse width and the number of frequeney steps. 
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Figure 20. Pulse Compression implementation. 


When the frequeney has been stepped through the whole pulse, the drops 
are moved in aeeordanee to the ehosen PRT. The proeess then starts the frequeney and 
range stepping again, eolleeting returns from the next pulse. The total power from eaeh 
pulse is summed over all frequeneies and then averaged for the number of range bins 
used. 


e. Plot of Phasors for Phase Change Comparison 

To visually illustrate the phase ehange for pulse pairs over time, 30 pulse 

pairs are presented in polar plots as phasors (see example Figure 21). The eurrent version 

of the simulator ean only plot 30 pulse pairs but uses all simulation pulse pairs to 
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calculate the average radial veloeity estimate. The estimated radial veloeities retrieved 
from eaeh pulse pair are plotted in a Histogram with a Gaussian fit for comparison, as 
presented in Figure 22. A three frequeney stepped scheme simulation (pg. 44) produced 
an average radial veloeity estimate of v)=-7.38 m/s and a standard deviation of 
(7^, = 4.68 m/s for a single run of 51 pulses calculated using Equations (2.32), (2.33) and 
(2.34) (see Figure 22). 


1^. v= -15.6291 m/s 2""^. v^= -12.6913 m/s 3'''^. v= -4.566 m/s 



Figure 21. The three first pulse pairs plotted in a Phaser plot. The input 
radial velocity was -10 m/s with a velocity spread of 4 m/s. 


Radial Velocity Histogram with Gaussian Fit 



Radial velocity [m/s] 


Figure 22. Radial Velocity Histogram with Gaussian Fit for 51 pulse 
pairs. Input velocity was -10 m/s with a velocity spread of 4 m/s. 
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f. Adding Spread to the Doppler Spectrum 

As described in a previous section, a spread of the velocity spectrum needs 
to be implemented to make differences in correlation possible. Applying a Gaussian 
random component to the wind velocity to introduce a controllable value of velocity 
spread was tested. This enables verification of the theory of the spread contributing to 
de-correlation. Results from running simulations, again using a PRF of 2000 Hz and no 
fall angle, for different spreads applied can be observed in Figure 23. The spread applied 
was, from left to right, <7„ = 4 m/s, = 3 m/s, and <7^, = 2 m/s, and the average 

simulator velocity spread estimates were = 4.05 m/s, <7^, = 2.94 m/s, and 
(7^, = 2.00 m/s again calculated using Equations (2.32), (2.33) and (2.34). 



Doppler FiequBncy [Hz] Doppler Frequency [Hz] Doppler Frequency [Hz] 

Figure 23. Test of Doppler spread, (a) input spread was 4 m/s, output 
4.05 m/s. (b) input spread was 3 m/s, output 2.94 m/s. (c) input spread was 2 m/s, 

output 2.00 m/s. 

The best de-correlation results were accomplished when the velocity 
spread was applied in two steps; first, a spread representing the spread of wind velocity 
over time, second, a spread representing the spread due to differences in drop sizes and 
other effects that will make the velocity of each drop vary over time. The fact that drops 
fall out of the resolution volume and new drops join adds an uncontrollable spread 
mechanism that can be examined by running several consecutive runs with the same 
input data. As the number of samples increases, the spread of the average velocity is 
expected to become smaller [3], which can be used for comparison. 
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D. SIMULATOR DESCRIPTION 


The simulator was built in modules developed to enable desired funetionalities. 
Initially the modules were separate funetions that a main program ealled for exeeution, 
but MATLAB works faster when only one long seript is used rather than funetions. The 
inputs are separated into Radar data and Rain parameters and are speeified in previous 
seetions. Part I of the simulator uses input Radar data to generate a Beam Resolution eell 
approximation and a Cube with an added margin to the Beam Resolution eell to allow for 
the rain drops to fall the through beam over time. Part II uses the Rain Parameters to 
generate a fundamental T Matrix eonsisting of: 

• Diameter Veetor and Number of drops in eaeh diameter interval using the 
Marshall-Palmer exponential approximation for drop size distribution. 

• The baekseatter RCSs for eaeh Diameter size 

-5 

• The total number of drops of eaeh diameter per m 

When the frequeney stepped version is used the RCSs are plaeed in a separate 
matrix, sinee RCS is frequeney dependent. In Part III all drops a plaeed randomly in the 
Cube developed in Part I and in Part IV the eoherent Eleetrie field return is ealeulated for 
every pulse. Between pulses all drops are moved in aeeordanee with their veloeities 
where veloeity has a random eomponent determined by . In Part V the final 
ealeulations are completed and the results are plotted. 


I. Part I 

To create the radar resolution cell based on the Radar data the volume of interest 
is retrieved by calculating the volume of the cone segment that represents the resolution 
cell (see Figure 9). Assuming an elliptical beam, the resolution volume can be 
approximated as [5] 






J 


(3.8) 


To allow for the water drops to fall through the resolution cell over all pulses, a 
margin is added to the resolution cell creating a box where all drops are initially placed. 
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Both volumes are defined in spherieal eoordinates to allow for plaeement, plots, and later 
eomparison. Also the volumes are registered. 


2. Part II 

This part produees a matrix of all drops that eontribute to the total refieetivity of 
the resolution eell. Using the diameter limits defined in Rain parameters, a drop size 
veetor is ereated. The limits ean be ehosen to ereate results elose to the integral form by 
pushing the upper limit of diameter to 15 mm (see results Figure 12), or to resemble more 
realistie values setting the max diameter to around 6 mm. Depending on what limits are 
ehosen, the aetual Z must be eomputed using the appropriate model, i.e. chose between 
Equation (2.30) and (2.31). The simulator allows for differences in resolution for 
different diameter ranges depending on what one wants to examine. The diameter vector 
is then used to create an RCS vector applying Equation (2.40), and a Drop Size 
Distribution (DSD) vector applying Equation (2.35) (Marshall-Palmer). To improve the 
speed of the simulator the number of drops is scaled down by a factor {Norm), which is 
used later to rescale the output values. The number of drops is rounded off to nearest 
integer using the MATLAB function ROUND. Einally a matrix is created to store all 
values computed in Part II 


matrix T round 


D, 


D, 


N(D,) ... N(D„) 


(3.9) 


and errors due to rounding off are calculated. The errors accounted for are RCS, Z, 
Number of Drops, and Rain Rate. In the stepped frequency and pulse compression 
versions the RCSs are stored in a separate matrix 



• 

■ ■ ^ 

res = 
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j 


(3.10) 


3. Part III 

The drops created in Part II are placed inside the Box created in Part I. Eor each 

diameter every drop is randomly placed inside the Box using a uniform distribution 
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function and its position is stored in spherical coordinates in a position matrix saved in a 
.mat-file 


( 


matrixpos = 


0 , 


yangle^ 


where angle is the angle off boresight. 


(PN 

0 , 

angle 


(3.11) 


4. Part IV 

In Part IV the Coherent Electric Field return is calculated for every pulse 
specified in the Radar Data. For every diameter, drops that are positioned inside the 
Resolution Cell are selected and based on matrices (3.9) and (3.11) the weighted electric 
field return is calculated using Equation (3.5) and summed for all drop diameters. Also 
the rain rate, for the specific Resolution Cell is calculated using Equation (2.39) 
approximated by a summation over all drop sizes. After the calculations, all drops are 
moved to a new position by applying the wind vector and terminal fall velocities 
specified in Rain parameters using Equation (2.37). All drops are assumed to move with 
the velocity of the wind vector with a Gaussian spread applied. The spread can be 
applied equally to all drops, so all drops move with the same velocity between pulses and 
the spread is applied over all pulses, or applied per diameter size so equally sized drops 
move at the same velocity between two specific pulses and the spread is applied over 
sizes. A third possibility is to apply the velocity spread on all drops individually creating 
a true spread although not fully realistic. The new drop positions are stored in matrixpos. 

In the stepped frequency version, two different PRTs are used together with a 
delay described in previous sections. After the first pulse has been registered in terms of 
I and Q return, using frequency I, the drops are moved according to the Doppler PRT 
{PRT_1). The phase is then registered completing the pulse pair before moving the drops 
again, this time using the Doppler PRT and the delay marking the change of frequency. 
Frequency 2 and frequency 3 are then applied in the same way as frequency I. When 
changing back to frequency I the simulator uses the second PRT {PRT_2) and subtracts 
the time elapsed due to the other frequencies and their transmitted pulse pairs. The return 
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values are stored in a matrix making sure the errors and the final estimates aeeount for the 
differenees in wavelengths and range between the returns. 

The pulse compression version steps through the number of stretched pulses 
specified in the Radar data. For every stretched pulse the chirp is approximated by 
discrete frequency steps using the frequency step size and the number of steps defined. 
Every frequency is associated with a specific range bin defined by the Resolution Cell 
depth and the number of frequency steps. After the return values for all frequencies and 
range bins have been collected, all drops are moved and a new pulse initiated. The time 
between pulses is defined by the PRT. 

5. Part V 

The final part of the simulation produces the weather signal parameter estimates 
and plots the results. The estimated parameters are 

• Average Power estimate, P , retrieved using the weighted electric field 
returns. 

• Reflectivity estimate, Z, retrieved from the average power estimate using 
(2.28) and compared with Z based on the simulator inputs using (2.27), 
evaluated by a summation and Equation (2.30), evaluated using the 
parameters Aq and A from the Marshall-Palmer DSD. 

• Average doppler frequency estimate, /^, retrieved using a pulse pair 
algorithm that extracts the phase-change over time (see Equation (2.32)). 

• Doppler spectrum width estimate, , retrieved taking the standard deviation 
of the doppler frequency data (see Equation (2.33) and (2.34)). 

Here the errors and the reduction factor are put back in the calculations. As a 
comparison, the ground truth for the respective parameter is derived using the integral 
form for calculation and plotted in the same graphs for clarity. Examples of output plots 
are presented in Eigure 24 where a simulation using three frequency steps was evaluated. 
The input parameters include frequency steps from 3 GHz to 3.002 GHz with a 1 MHz 
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frequency step, a Doppler PRF of 5,000 Hz, and the second PRF of 300 Hz. Further, the 
radial velocity was -10 m/s with a velocity spread of 4 m/s, which should produce a 


p ~ 0.246 , using 7 =- 

"300 


s. 


The average Doppler frequency estimate =-150 Hz 


(Figure 24 (a)) is about 25 % from the expected Doppler frequency of = -200 Hz. 

Only 51 sample pairs were used to estimate the average Doppler, which could explain the 
difference between the estimate and the true value. The simulation is a point estimate, 
which means that running consecutive trials would better show the behavior of the 
velocity estimate. The offset can also be the effects of adding new drops and loosing 
others in the resolution cell as explained in previous sections. The plot of Estimation of 
Z, (c) presents results close to the outcome expected using the mathematical formulas for 
comparison. The correlation value of 0.246 can explain why the value of the estimate Z 
is so close to the integral form, as de-correlated samples will yield an estimate with a 
smaller variance. The example output plot of correlation coefficient, (d), although close 
to the expected value from the above described simulation, is from another simulation 
running only one PRF. The input parameters were: 


• Wavelength, A 0.1 m 


• PRF 300 Hz 


• Wind -10 m/s in x-direction (=radial velocity) 

• Velocity spread 4 m/s 

The expected correlation, using Equation (2.15), is ~ 0.246 , which is very close 
to the simulator output value. 
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(a) 


(b) 


Doppler frequency ftjr all pulse pairs for PRF= 5000 Hz 



Power Return for all 51 pulses. Frequency stepped from 3 GHz to 3.002 GHz with 1 MHz steps. 



(c) 


Estimated Zper pulse 



Return Voltage Correlation 



Figure 24. Example output plots from simulator, (a) Doppler, (b) power, 
(c) reflectivity, and (d) correlation coefficient. Note that (d) is from a different 

simulation than (a) - (c). 
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IV. ANALYSIS OF APPLICATIONS AND DISCUSSION OF 

RESULTS 


A. IS IT POSSIBLE TO BUILD A WEATHER RADAR SIGNAL 

SIMULATOR BASED ON A PHYSICAL MODEL OF A SPATIAL 

REGION CONTAINING RAINDROPS? 

The physical model built is based on the Rayleigh scattering approximation and 
the number of scatterers and their sizes are specified using the Marshal-Palmer drop size 
distribution. The drops are approximated by spheres and are uniformly distributed in a 
volume defined by the radar data and the return signal is calculated by summing the 
electric field. The model accounts for varying terminal velocities for different drops sizes 
and allow input for wind, which is represented by a 3-dimesional vector applied to every 
drop. The wind vector is implemented with a Gaussian spread among drops enabling 
investigation of correlation effects, which can be applied within or between pulses. 

The results from running the simulator show that using the model developed will 
generate data following expected theoretical distributions. The simulator treats the rain 
drops as point targets adding the electric fields to generate the RCS. The results verify 
the assumption that the average return from rain can be approximated by adding the 
RCSs of all drops. The simulator can be adapted to explore effects of changes in DSD. 
Since the simulator distributes the drops randomly within the whole volume, including 
the resolution volume and its margins, the average distribution for sub-cells does not fully 
follow that of Marshal-Palmer. The effects of this should be further explored. 

The current version of the simulator does not take into account attenuation, 
mutual coupling, and multiple scattering. Neither losses nor noise are introduced and 
dual polarization is left for future work. Although simplified, the weather radar signal 
simulator based on a physical model of a spatial region containing raindrops produced 
from this research work will allow for further development improving realism. 
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B. CAN REALISTIC MOTION THAT DEPENDS UPON THE TYPE OF 

WEATHER SYSTEM BE IMPARTED TO THE RAINDROPS? 

This work has prioritized the eleetromagnetie and radar aspeets and the 
meteorology portions have been kept to a minimum. The simulator version presented 
uses a veetor model to introduee motion to the raindrops. This ean be adapted to impart 
realistie motion. To do that, realistie models need to be aequired in veetor form allowing 
implementation in the simulator. The eurrent version of the simulator measures only one 
resolution eell exploring effeets of de-eorrelation using different schemes. However, the 
simulator can be extended to cover a larger volume including several range bins. Using 
velocity spread as an input parameter makes it possible to simulate spectrum widening 
effects without having to create a detailed physical model of each spreading factor. 
There is no doubt that using a vector model would make it possible to implement realistic 
motion to the raindrops although this has not been done in this version of the simulator. It 
was beyond the scope of the thesis to fully implement realistic motion models. 


C. CAN THE OUTPUT OF THE WEATHER RADAR SIMULATOR BE 
USED TO STUDY SPATIAL AND TEMPORAL SAMPLING SCHEMES 
FOR CONSTANT OR STEPPED FREQUENCY SAMPLING PULSES? 

Two methods have been implemented focusing on temporal sampling; namely 

adapting Pulse Repetition Frequency so the drops have time to reshuffle between pulses, 

and Frequency shifting between sample pairs. A simulation was run with the following 

data (for details see appendix C): 


• 

Wavelength, X 

0.1 m 

• 

PRF 

300 Hz 

• 

Wind 

-5 m/s in x-direction (=radial velocity) 

• 

Velocity spread 

2 m/s 

To 

determine a suitable PRF the sample time threshold provided in Equation 


(2.18) was used. The average Z of the simulator data lies within 0.5 dB of the integral 
form which is also true for the value retrieved using the power samples (see Figure 25). 
The reflectivity estimate based on sample power lies approximately 1 dB below the 
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estimated average refleetivity, Z and 0.5 dB below the ground truth, whieh ean be 
explained by the approximation of the radar beam where only drops inside the half power 
beamwidth are eonsidered. The drop position within the radar beam will also affeet the 
gain level, whieh makes the position of the drops vital. A eonstant gain approximation 
might result in a more aeeurate approximation. The ehosen PRF led to a level of 
eorrelation of p ~ 0.543, whieh is just below the threshold of 0.6 but still enough to 
make it diffieult to estimate radial veloeity with reasonable aeeuraey. The sample pairs 
are not strongly eorrelated whieh makes an estimate possible but the varianee will eause 
errors partly due to the laek of eorrelation, and partly beeause there are too few samples. 

A new simulation was run this time using a rapid seanning seheme with stepped 
frequeneies as deseribed in a previous seetion. The new simulation was run using the 
same inputs but adding another PRF of 5000 Hz to improve the Doppler estimation 
eapability. 

The result were very similar to the first run using only time to aequire de- 
eorrelation, however, eomparing the total time to produee 50 samples differ signifieantly, 
still aeeomplishing a Z-value elose to the integral form (Figure 26). The reason for this is 
that the 300 Hz PRF only needs to be used to separate equal frequeneies while stepping 
between frequeneies ean be done very quiekly. Three de-eorrelated samples ean now be 
eolleeted in almost the same amount of time as previously only one pulse eould be. The 
pulse pair seheme also enables an estimation of radial veloeity seen in Figure 27. The 
spread was larger than the input of 2 m/s whieh ean be explained by the faet that the 
spread was applied both between pulses and within pulses between different drop sizes. 
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Estimated Z per pulse 
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Figure 25. Temporal Sampling, PRF variation. PRF 300 Hz. Velocity 
Spread 2 m/s between pulses and within pulses for drops of different sizes. 


The output from the latest simulation presented an average radial veloeity of 
-5.25 m/s and a veloeity spread of 4.6 m/s (using Equations (2.33) and (2.34)) whieh is 
elose to the expeeted values. 
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Estimated Z per pulse 
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Figure 26. Temporal Sampling, Stepped Frequency. PRFi 5000 Hz, PRF 2 
300 Hz. Velocity Spread 2 m/s between pulses and within pulses for drops of 

different sizes. 
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Figure 27. Velocity histogram with Gaussian Fit, Temporal Sampling 
using stepped frequency. The input velocity was -5 m/s. (J, = 4.6 is estimated 

using Equations (2.33) and (2.34). 
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To analyze the simulator behavior over time for several eonsecutive runs, a 
simulation series was performed using 50 runs and the results were analyzed for different 
number of pulses. In [2] the intensity probability distribution for the averaged return for 
different numbers of pulses was plotted as in Figure 28, where the dependenee of the 
number of pulses can be seen. 



Figure 28. Probability distribution of Jk, which is the return intensity. 
The intensity average is (From [2]). 

Running the simulator using a A = 0.1 m, and a PRF of 200 Hz, which represents 
a correlation coefficient of about 0.04 (Equation (2.15)) considering a velocity spread of 
4 m/s, renders results seen in Figure 29. The figures resemble each other in shape but for 
A/^ = 100 samples, the peak for the simulator distribution is not as high as the result in 
Figure 28. As the number of pulses goes up the distribution become more Gaussian, 
which is expected for any distribution if the number of samples increases. Also, the 
spread decreases as the number of pulses increases. Having completely independent 

samples will change the average velocity spread as <7- = , where k is the number of 
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pulses or samples. The results are somewhat expeeted sinee the distribution of the power 
return is expeeted to be exponential [3]. The standard deviation of an exponential 
distribution is expeeted to be equal to the mean of the same, whieh the results verily. 
However, if samples are independent the spread should deerease with a faetor of 1 / 4k 

as k inereases. For 2 pulses the expeeted value equals 1/ V2 = 0.71, whieh matehes the 
simulator result. For 10 pulses the value is 0.316, while the simulator shows 0.41. 
Finally for 100 pulses the value is 0.10. The simulator shows 0.21, which is more than 
the double what is expected. This suggests that not all samples can be considered 
independent. Why this occurs when a sufficient spread is applied cannot be explained at 
this point but is of relevance for future investigations. 

Histogram 



Figure 29. Histogram showing distribution of intensity (power) from 200 
consecutive runs. Number of pulses averaged is 1, 2,10 and 100. p ~ 0.04 so 

independent samples are expected. The relative standard deviation (/ P) in 
each case is, 1 pulse=0.99,2 pulses=0.74,10 pulses=0.41, and 100 pulses=0.21. 
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D. CAN THE WEATHER RADAR SIMULATOR BE USED TO STUDY THE 
USE OF PULSE COMPRESSION AND RANGE AVERAGING AS A 
MEANS OF RAPIDLY OBTAINING INDEPENDENT SAMPLES? 

A simulation was performed using the simulator version adapted to pulse 

eompression as described in a previous section. The resolution volume differs from the 

previous simulations because of the stretched pulse making the range dimension of the 

resolution cell 10 times larger. All other data was set as before, i.e. PRF of 300 Hz, wind 

speed of 5 m/s and a velocity spread of 2 m/s. The number of pulses is still 50 but they 

are used differently letting every pulse cover a range bin equal to previous runs. The new 

stretched pulse consists of 10 sub-pulses of different frequencies creating the chirp. Each 

new pulse is summed over all frequencies stepped over the whole resolution cell and then 

averaged. Effectively only 5 pulses are used, but since each pulse consists of 10 

frequency steps, the average value is based on 50 samples. 


Estimated Z per pulse 



Figure 30. Pulse Compression simulation. 5 pulses, Frequency stepped 
from 3 GHz to 3.0018 GHz with 0.2 MHz steps for each pulse. 


The estimate of Z based on averaging lies 1.0 dB above the result using the 
integral form, which indicates good accuracy. The smaller deviation between samples 
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compared with previous simulations ean be explained by the faet that eaeh value really 
represents an average of 10 range bins. The higher level ean be explained by the method 
used to average. The volume of the different resolution bins inereases moving away from 
the radar so if this is not taken into aeeount when averaging, the values will end up 
deviated from the expeeted value. Another explanation is the faet that the spread in 
frequeney is not enough to ensure suffieient de-eorrelation to get independent samples. 
Running the simulator for 8 stretched pulses generated a result very elose to the integral 
form (0.4 dB) something that needs further investigation for verifieation. 

The results show the possibility to use the simulator to study the use of pulse 
eompression and range averaging as a means of rapidly obtaining independent samples. 

E. CAN THE OUTPUT OF THE WEATHER RADAR SIGNAL SIMULATOR 
BE USED TO PRODUCE DISPLAYS OF PEDAGOGICAL INTEREST? 

During the progress of building the simulator, several displays have been 

developed and used to assist with the evaluation of eaeh module. Onee the simulator 

output data has been produeed it is all numbers making it possible to post proeess in 

many different ways, displays being one. To be able to use the simulator in teaching, 

visual aids presenting the signal pulse-to-pulse were developed as well a phasor 

representation displaying the phase ehange over time between pulse pairs. Adding a 

Graphieal User Interfaee, GUI, would probably enhanee the utility the simulator as a 

pedagogieal tool for investigating the effeets of different weather phenomena and 

variation of input values for both radar and weather parameters. Nevertheless, the output 

of the existing simulator ean be used to produee displays of pedagogieal interest. 

F. CAN THE OUTPUT OF THE WEATHER RADAR SIGNAL SIMULATOR 
BE USED TO STUDY THE PERFORMANCE OF ESTIMATORS FOR 
THE FIRST THREE DOPPLER MOMENTS? 

The first versions of the simulator only allowed for ehanges in a single PRF 
making it possible to vary time between pulses effeeting de-eorrelations. The simple 
form of return signal produeed samples equally spaeed in time making it simple to 
perform post proeessing in the form of auto eorrelation and FFT. The speetrum was 
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plotted and it was possible to derive an estimated radial veloeity of the drops from the 
average Doppler frequeney. When the eomplexity of the return signal inereased and the 
simulator allowed for multiple PRFs and delays, the signal needed manipulation to enable 
FFT or autocorrelation. Initial trials padding the signal with zeros making the samples 
equally spaced in time were performed but adding time to a signal demands more signal 
processing to remove effects of folding and aliases. Digital signal processing is beyond 
the original scope of the thesis which is why this interesting trail was abandoned. There 
are no test runs verifying that the simulator can be used to study the performance of 
estimators for the first three Doppler moments but that does not imply that the simulator 
cannot be used for this application. Having a complex return signal allows for post 
processing of many sorts, as mentioned earlier. 


G. CAN THE OUTPUT OF THE WEATHER RADAR SIGNAL SIMULATOR 
BE USED TO STUDY THE UTILITY OF ESTIMATORS FOR HIGHER 
ORDER DOPPLER MOMENTS? 

Since the estimation of the lower order of estimators was not performed, the 
higher order estimators have not been investigated. 
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V. CONCLUSIONS AND RECOMMENDATIONS FOR FUTURE 

DEVELOPMENT AND WORK 

A. CONCLUSIONS 

The objective of the research was to build a weather radar signal simulator 
delivering I and Q channel outputs from a physical representation of rain drops. A 
previously developed simulator was initially evaluated and verified serving as the basic 
building block when developing the final version. The Radiation Integrals were used to 
derive the electromagnetic scattering from rain drops approximated as dielectric spheres. 

Without accounting for multiple scattering, attenuation, or coupling effects, the 
scattered electric field from each drop is summed for every pulse providing the total RCS 
of the rain in the radar resolution cell as well as a complex voltage return representation. 
Using the Marshall-Palmer drop size distribution, an initial position for every drop is 
determined within an extended volume enabling the drops to fall through the radar beam 
for all samples. Between pulses the drops are moved according to their respective 
terminal velocities based on the Atlas-Ulbrich approximation and wind. To allow for 
width in the velocity spectrum, a Gaussian velocity spread is applied to the drops 
between and within pulses. The simulator power return is weighted using the Radar 
Range Equation for point targets, considering all drops as point targets, and evaluated as 
a volume target using the Weather Radar Range Equation. 

Several different functionalities have been implemented allowing for stepped 
frequencies, multiple PREs, pulse compression using a chirp, and variation of radar and 
weather input parameters. Post processing capabilities include autocorrelation and PET 
(only for single PRE version); computation of weather parameter estimates such as 
reflectivity factor, Z; average doppler, radial velocity, and velocity spread; pedagogical 
plots including a Phaser plot of phase change over time, a velocity histogram, power and 
reflectivity for each pulse over time. 

This research has shown that building a weather radar signal simulator based on a 
physical model of a spatial region containing raindrops is possible. The benefit of using a 
physical model is improved realism and the ability to study effects of scientific interest. 
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Also, the model permits isolation of effeets making it possible to explore the results when 
varying some parameters while keeping others fixed. Using a veetor model to input wind 
and other effeets on the rain drops makes it possible to impart realistie motion to the 
drops although this has not been done fully in this version of the simulator. This was 
beyond the seope of this thesis but is a reeommended area for further study. Evaluating 
the result from averaging radial veloeity rendered an interesting finding namely the 
differenee between using EFT, autoeorrelation, and pulse pair proeessing to estimate 
radial veloeity. Using the EFT or autoeorrelation skewed the result towards the bigger 
drops beeause the amplitude information is kept in the proeess. On the other hand with 
the pulse pair algorithm, only the phase information is used, resulting in a “true” average 
veloeity for all drops regardless of size. Intuitively smaller drops would be more 
suseeptive to wind than bigger drops which would emphasize the need for pulse pair 
algorithms over the EFT or autocorrelation for better wind velocity estimates. No 
theoretical investigation or look for research in this area has been performed but is 
suggested for future studies. Different implementations of velocity spread have been 
tested leading to variations in results. The final version of the simulator allows for spread 
both between pulses, which represents how the average velocity is spread between 
samples; and within pulses with a size dependent spread, which represents the spread 
within the resolution cell. 

The title of this research indicates where the priority lies, namely developing the 
simulator to examine sampling rates and scanning schemes. The results show that it is 
possible to use the weather radar simulator to study spatial and temporal sampling 
schemes for constant or stepped frequency sampling pulses. The research aimed at 
building the radar and not to fully explore the effects of different rate and schemes, but 
the results from testing and verification indicate that using methods to more rapidly 
obtain independent samples decreases the time to obtain accurate estimates of weather 
signal parameters. By introducing three frequency steps, the time to acquire accurate 
data was reduced by a third, still being able to acquire Doppler information as well as 
reflectivity. However, more simulations need to be performed comparing results when 
parameters are varied to fully explore the benefits of frequency stepping and its 
limitations. 
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By implementing a pulse eompression eapability, using frequeney stepping and 
range averaging, results were provided showing the ability to rapidly obtain independent 
samples using pulse eompression. Inereasing the radar resolution eell slows down the 
simulator, whieh ean be eompensated by inereasing the A^orm-faetor (see previous 
ehapters for details). Again further testing is suggested exploring the effeets of varying 
parameters. 

Post proeessing the signal samples provides wide-ranging opportunities for 
different displays of pedagogieal interest. Most plots have been developed to eompare 
and verily results but some plots have been developed only to display pedagogieally 
interesting phenomena. Being able to follow phase ehange over time is an example of the 
latter. Implementing a GUI will probably be the greatest improvement in this area. 

Having an aeeurate signal allows for more than just plots of results, it opens up 
possibilities for post proeessing sehemes beyond the scope of this work. Previously 
indicated processing includes filtering, FFT and autocorrelation. In the introduction, 
cross scientific benefits were mentioned signifying the use of knowledge from one area in 
another. Having a weather signal consisting of rain echoes can be used to show or 
investigate filtering to remove weather clutter. 

The initial intentions included ideas of improving the simulator speed by 
optimizing algorithms and scripts. This has not been accomplished. However, the 
simulator runs a simulation locally on a PC with an Intel Pentium 4 processor 550 
(3.4GHz, IM, 800MHz FSB, NTFS Memory 4.0GB DDR2 Non-ECC SDRAM, 400Hz 
(4DIMM)), in just a few minutes, depending on complexity and input parameters. 

B. RECOMMENDATIONS FOR FUTURE RESEARCH 

The simulator is ready to do its job, setting more than one path for future work. 
The focus in this research has been to develop a tool for research, which leaves the 
weather research parts for future work. Two main areas of development and research are 
proposed: 

• Examine sampling rates and scanning schemes for fast volumetric scanning 
and update 
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Simulator development 




1. Examine Sampling Rates and Scanning Schemes for Fast Volumetric 

Scanning and Update 

Weather Radar Signal Simulator should be run to examine weather signal 
parameter estimation sensitivity when varying radar parameters and preeipitation 
estimator models where changing the drop size distribution would be of interest. 
Comparison should be made to real life measures as well as theoretical models. To 
complement this work weather models should be explored to find the physics behind 
different weather phenomena enabling implementation in the simulator. Alterations in 
the simulator might be necessary to facilitate this implantation. 

2. Simulator Development 

Even though the final version of the simulator has several interesting 
functionalities many interesting effects and functionalities can be added. 

a. Adding a GUI 

The main priority developing the simulator should be to create a GUI making the 
tool more user-friendly. 

b. Modeling Non-Spherical Drops 

The current simulator approximates each drop, regardless of size, as a 
dielectric sphere. Research show that drops flatten [9] while falling so an oblate spheroid 
approximation would be more appropriate. Work done by Bringi and Chandrasekar [14] 
provide good mathematical tools enabling implementation in the simulator. 


c. Adding Dual Polarization 

Implementing a spheroid approximation makes dual polarization 
measurements possible and new estimates can be computed such as differential 
reflectivity. The complexity of this can be overwhelming considering effects like cross 
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talk and multiple scattering. This is a problem for scientists or students who are 
searching for a real challenge. 

d. Post Processing for Higher Order Moments 

The output signal can be post processed to produce the power spectrum, 
which was done in the early versions that had uniform time spacing between samples. 
Applying digital signal processing to the output signal, retrieving the power spectrum 
would make it possible to explore higher order moments of the spectrum and to evaluate 
their utility as estimators. 

e. Linking Simulator Input to Storm Physics 

As mentioned, storm models should be explored to determine the physics 
associated with different weather phenomena. The physics can be used to determine the 
wind velocity and velocity spread that should be used as simulator input parameters. 
This will add to the realism of the simulator and will make it possible to explore effects 
on the radar signal caused by different weather phenomenon. 

f Adding Noise to the Weather Signal 

Real weather signals will always be contaminated by some noise. Noise 
from all sources should therefore be added as a set of inputs. The benefit of this would 
be the ability to study the degradation caused by noise. 


g. Adding Mutual Coupling 

Mutual coupling was examined but needs more attention and different 
techniques to enable implementation into the simulator. For large rainfall rates, raindrops 
may be spaced at distances that result in mutual coupling. The total RCS in this case will 
be less than that computed assuming no coupling. This effect needs further investigation 
to determine what storm types and rain rates result in coupling that renders the 
independents scattering assumption invalid. 
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APPENDIX A. DERIVATION OE RCS FOR SINGLE DROP 


A water drop can be approximated as a small dielectric sphere. Let the sphere 
radius be a, and the permittivity f, and permeability , also let the sphere be centered at 
the origin. An incident plane wave is present 

Ei = zE,e-^^ 

The electric field at a point in space can be expressed as 

E = Ei + Es 

Using the volume equivalent principal, the current J can be written as 

1 = jco{e-e^)E 

Letting a»/I, E is then approximately E = zE^ (constant inside). Using 

Es = L'l can be derived as 
3(0£q 


zEi 


E — Ei + Es — Ei + 


jJ 

Zcoe^ 


zEi = Ei + j 


ico{e^e,-e^)E 


Zcoe. 


3 


zE\ = 



3Ei 

(f.+2) 


Now solving for J 
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3Ei 
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Finally putting the current in the radiation integral 


E{r,0,(/)) = e Je^'‘^dv',g = 0 


47rr 
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Using this in the formula for RCS renders, (letting a — D/2 ) 


cr^ = 4;rr^ 
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APPENDIX B. MICROWAVE STUDIO SIMULATION SETUP 


First the target(s) was created. All targets were spheres of 4 mm in diameter with 
the specifications of a water drop. The incident ray was a plane wave in the z-direction 
with polarization in the x-direction or in the y-direction. When adding the second drop, 
the drop was placed side by side with the first drop, i.e. in the x-direction. The plots are 
farlield plots in spherical coordinates. The mesh resolution was set to 20 lines per 
wavelength. 



Material = waterdrop 

Type B Lossy metal 

hue = 1 

El. cond. = 0.0001 [5/m] 

Rho = 1 rkg/m'‘3] 


Figure 1. Target. Water drop. 
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Figure 2. Incident Plane for the one-drop simulation 



Figure 3. Type of plot and definition of angles. 
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Figure 4. Resolution of the Meshed grid. 
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Figure 7. Incident Plane for the two-drop simulation 
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Figure 8. Test of noise by adding extra volume when calculate RCS of one drop. 
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Figure 9. Plot of results adding extra volume when calculating RCS for one drop in the 0 - 

plane for ^ = 0. 
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APPENDIX C. 


Radar Data 

P_t=10e3; 

G_dB=20; 

G=10^(G_dB/10); 

angle_3db=0.20; 

pulsewidth=0.20; 

elevation=3; 

margine=5; 

azimuth=0; 

number_pulses=5 0; 

PRT=l/300; 

distance=0.5; 

wavelength=0.1; 

k=2*pi/wavelength; 

Rain parameters 

diam_low=0.05; 

diam_lim=l; 

diam_mid=5; 

diam_high=15; 

resolutionl=0.08; 

resolution2=0.3; 

resolution3=0.8; 

sort_approx=l; 

kw=0.93; 

rain_rate=50; 

rain_anglel=0; 

rain_angle2=0; 

wind_vector=[-5,0,0]; 

spread=2; 

Norm=100; 


SIMULATION SETUP EOR THE WEATHER 
RADAR SIMULATOR 


% Transmitted power 
% Maximum gain in dB 
% Maximum gain 

% 3 dB Beam angle [degrees], (half half power beamwidth) 
% Pulsewidth [mieroseeond] 

% Elevation angle [degrees] 

% margin for resolution oell[m] 

% Azimuth angle [rad] 

% Number of pulses 
% Pulse Repetition Time [seeonds] 

% Range from the radar to the resolution eell [Km] 

% Wavelength [m] 


% Lowest diameter of Preeipitation [mm] 

% Divider between low and mid [mm] 

% Divider between mid and high [mm] 

% Highest diameter of Preeipitation [mm] 

% Lower diameter resolution [mm] 

% Higher diameter resolution [mm] 

% Highest diameter resolution [mm] 

% RCS approximation model; [l]=Rayleigh 
% abs(K)^2 approx for water 
% Rain Rate [mm/h] 

% degrees Rain Angle (Phi) 

% degrees Rain Angle (Theta) 

% Wind speed veetor (x,y,z) [m/s] 

% Wind speed spread [m/s], (all direetions) 
% Normalization eoeffieient 
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APPENDIX D. CODE DEVELOPED IN MATLAB. 


There are eurrently four versions of the simulator: 


1) The basie version, WR basie, that will run with one PRF. The plots inelude 

• Doppler per pulse pair with average Doppler, 

• Radial veloeity histogram with Gaussian fit 

• Return power with average Return power 

• Refleetivity estimations with Ground Truth for eomparison 

• Voltage Correlation funetion with estimation of p{t = T^) 

• Power speetrum of return signal 

2) Dual PRF version, WR_dual_prf, that will run with two PRFs. The plots 

inelude 

• Doppler per pulse pair with average Doppler, 

• Radial veloeity histogram with Gaussian fit 

• Return power with average Return power 

• Refleetivity estimations with Ground Truth for eomparison 

• Phasor plots for the 30 first pulse pairs. 

3) Frequeney step version, WR_fq_step, that will step through several frequeneies 
and transmit a pulse pair on every frequeney. One PRF eontrols the time between pulses 
in the pulse pair and one PRF eontrols the time between equal frequeneies. The plots 
inelude 

• Doppler per pulse pair with average Doppler, 
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• Radial velocity histogram with Gaussian fit 

• Return power with average Return power 

• Reflectivity estimations with Ground Truth for comparison 

• Phasor plots for the 30 first pulse pairs. 


4) Pulse compression version, WR_pulse_com, that divide the resolution volume into 
range bins and averages over bins and pulses. The chirped pulse is approximated by 
discrete stepped frequencies. The plots include 

• Return power with average Return power 

• Reflectivity estimations with Ground Truth for comparison 

1) WR_basio 


% Weather Radar Signal Simulator, developed version 
% Ulf Sehroder 
% July, 2005 

% - 


% Main program for The NPS/Marielle Gosset Weather Radar Signal Simulator. 


% - 

% DESCRIPTION 
% Part I 

% From Radar Data: 

% - Generate Beam Resolution eell approximation 

% - Generate Cube adding a margine to the Beam Resolution eell to allow for 
% the rain drops to fall trough the through beam over time. 

% 

% Part II 

% From Rain Parameters: 

% Generating the fundamental T Matrix eoneisting of 
% - Diameter Veetor and Number of drops in eaeh diameter interval 
% using the Marshall-Palmer exponential approximation. 

% - the baekseatter RCSs for eaeh Diameter size 
% - the total number of drops for eaeh size per m''3 

% Part III 

% Randomly plaee all drops in eube 
% Part IV 

% For every pulse the eoherent Eleetrie field return is ealeulatet. Between 
% pulses all drops are moved. 


% Part V 
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% Calculate and plot results 

% - 

clear all 
close all 
clc 

% Input Data 

c=3e8; % Speed of light [m/s] 

% Radar Data 

P_t= 10e3; % Transmitted power 

G_dB=20; % Maximum gain in dB 

G=10^(G_dB/10); % Maximum gain 

angle_3db=0.20; % 3 dB Beam angle [degrees], (half half power beamwidth) 

pulsewidth=0.20; % Pulsewidth [microsecond] 

elevation=3; % Elevation angle [degrees] 

margine=5; % margin for resolution cell[m] 

azimuth=0; % Azimuth angle [rad] 

number_pulses=40; % Number of pulses 

PRT=l/200; % Pulse Repetition Time [seconds] 

distance=0.5; % Range from the radar to the resolution cell [Km] 

wavelength=0.1; % Wavelength [m] 

k=2*pi/wavelength; 

% Rain parameters 

diam_low=0.05; % Lowest diameter of Precipitation [mm] 
diam_hm=l; % Divider between low and mid [mm] 

diam_mid=5; % Divider between mid and high [mm] 

diam_high=15; % Highest diameter of Precipitation [mm] 

resolutionl=0.08; % Lower diameter resolution [mm] 

resolution2=0.3; % Higher diameter resolution [mm] 
resolution3=0.8; % Highest diameter resolution [mm] 
sort_approx=l; % RCS approximation model; [l]=Rayleigh 
kw=0.93; % abs(K)^2 approx for water 

rain_rate=50; % Rain Rate [mm/h] 

rain_anglel=0; % degrees Rain Angle (Phi) 

rain_angle2=0; % degrees Rain Angle (Theta) 
wind_vector=[-10;0;0]; % Wind speed vector (x,y,z) [m/s] 
spread=l; % Wind velocity spread [m/s], (all directions) 

spread_2=4; % Drop velocity spread [m/s], 

Norm=100; % Normalization coefficient 

% Directional cosines 

u=sin(pi/2-elevation*pi/180)*cos(azimuth*pi/180); 

v=sin(pi/2-elevation*pi/180)*sin(azimuth*pi/180); 

w=cos(pi/2-elevation*pi/180); 

% - 

% Part I, Generate the resolution volume 

% - 

% Approximation of the resolution cell of the pulse circular box, 

% for any azimuth angle 

% - 

% INPUT 
% - distance 
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% - angle_3dB 
% - pulsewidth 
% - elevation 
% - margine 
% - azimuth 

% - 

% OUTPUT 

% - volume resolution: eomputed volume of the resolution eell 
% - volume box: volume of the box where the drops are put 
% - box: the 8 reetangular eoordinates of the box 
% - eoord vol res: the eoordinates of the resolution eell where the drops 
% are taken. 

% - 


% - 

% Creating the eireular beam resolution eell 
range_res=(e*pulsewidth* le-6)/2; 

% Range resolution (et/2) [m] 

% small angle approximation allows for following 
angle_rad=(angle_3 db *pi)/180; 

% 3 dB beamwidth [rad] 
radius 1=(tan(angle_rad)* (distanee *1000)); 

% Radius of beam at R=distanee 
radius2=(tan(angle_rad) *((distanee * 1000)+range_res)); ... 

% Radius of beam at R=distanee + range resolution 

surfaee 1 =pi *radius U2; 
surfaee2=pi *radius2^2; 

% volume of eone = base surfaee times height for all three 
volumeeone 1 =(surfaee 1 *(distanee* 1000))/3; 
volume_eone2=(surfaee2*((distanee*1000)+range_res))/3; 
volume_resolution=(volume_eone2-volume_eonel); % Beam Cone Volume 

% eomputation of the resolution eell eoordinates, spherieal eoordinates 
% elevation in radians 

elevation_rad=(elevation*pi)/180; % Elevation [rad] 

azimuth_rad=(azimuth*pi)/180; % Azimuth [rad] 

% Matrix used when eomparing wheather a drop is inside or outside 
% beamwidth 

eoord_vol_res=[distanee*1000,(distanee*1000)+range_res;elevation_rad-angle_rad,... 

elevation_rad+angle_rad;azimuth_rad-angle_rad,azimuth_rad+angle_rad]; 

% - 


% - 

% Creating the box 

% eomputation of the box and its eoordinates 

eoinl=[azimuth_rad+angle_rad,elevation_rad+angle_rad,distanee* 1000]; 

eoin2=[azimuth_rad-angle_rad,elevation_rad+angle_rad,distanee* 1000]; 

eoin3=[azimuth_rad+angle_rad,elevation_rad-angle_rad,distanee *1000]; 

eoin4=[azimuth_rad-angle_rad,elevation_rad-angle_rad,distanee*1000]; 

eoin5=[azimuth_rad+angle_rad,elevation_rad+angle_rad,distanee* 1000+range_res]; 

eoin6=[azimuth_rad-angle_rad,elevation_rad+angle_rad,distanee*1000+range_res]; 

eoin7=[azimuth_rad+angle_rad,elevation_rad-angle_rad,distanee*1000+range_res]; 
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coin8=[azimuth_rad-angle_rad,elevationrad-anglerad,distance* lOOO+rangeres]; 


[xl,yl,zl]=sph2cart(coinl(l),coinl(2),coinl(3)); 
[x2,y2,z2]=sph2cart(coin2(l),coin2(2),coin2(3)); 
[x3 ,y3 ,z3]=sph2cart(coin3 (l),coin3 (2),coin3 (3)); 
[x4,y4,z4]=sph2cart(coin4(l),coin4(2),coin4(3)); 
[x5,y5,z5]=sph2cart(coin5(l),coin5(2),coin5(3)); 
[x6,y6,z6]=sph2cart(coin6(l),coin6(2),coin6(3)); 
[x7,y7,z7]=sph2cart(coin7 (l),coin7 (2),coin7 (3)); 
[x8,y8,z8]=sph2cart(coin8(l),coin8(2),coin8(3)); 

xmax=max([xl,x2,x3,x4,x5,x6,x7,x8]); 

xmin=min([xl,x2,x3,x4,x5,x6,x7,x8]); 

ymax=max([yl,y2,y3,y4,y5,y6,y7,y8]); 

ymin=min([yl,y2,y3,y4,y5,y6,y7,y8]); 

zmax=max([zl,z2,z3,z4,z5,z6,z7,z8]); 

zmin=min([zl,z2,z3,z4,z5,z6,z7,z8]); 


% the box with its margin 

coin_l=[xmin-margine,ymin-margine,zmin-margine]; 

coin_2=[xmin-margine,ymax+margine,zmin-margine]; 

coin_3=[xmax+margine,ymax+margine,zmin-margine]; 

coin_4=[xmax+margine,ymin-margine,zmin-margine]; 

coin_5=[xmin-margine,ymin-margine,zmax+margine]; 

coin_6=[xmin-margine,ymax+margine,zmax+margine]; 

coin_7=[xmax+margine,ymax+margine,zmax+margine]; 

coin_8=[xmax+margine,ymin-margine,zmax+margine]; 


box=[coin r,coin_2',coin_3',coin_4',coin_5',coin_6',coin_7',coin_8']; 

% calculate the volume of the box 

volume_box=abs((xmax+margine-(xmin-margine))*(ymax+8-(ymin-margine))*... 

(zmax+margine-(zmin-margine))); 


cube=[coinr,coin_2',coin_3',coin_4',coinr,coin_5',coin_6',coin_2',... 

coin_6',coin_7',coin_3',coin_7',coin_8',coin_4',coin_8',coin_5']; 

% - 


% - 

% Part II, Generate the fundamental T Matrix 

% - 

% This part produces a matrix of all drops that contribute to the total 
% reflectivity of the resolution cell 

% - 

% First all sizes are created, with two vectors of different size resolution. 

% The size vectors represents all diameters. Next the backscatter RCS for each 
% size is computed. 

% Using the Marshall-Palmer exponential approximation, the number of drops of 
% each diameter is computed. 

% - 

% INPUT 
% - diam low 
% - diam lim 
% - diam high 
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% - resolution! 

% - resoIution2 
% - sort approx 
% - wavelength 
% - Kw 
% - rain_rate 

% - 

% OUTPUT 
% - number drops 
% Matrix T: 

% - first row all sizes 
% - seeond row all the baekseatter RCSs 
% - third row the total number of drops for eaeh size per m^3 
% the eolumns starting from the smallest size 

% - 

% Create dropsize veetor 

% Preeipitaion Diameter vector (lowest<D<divider) 
vector_size_ 1 =diam_low resolution 1 :diam_lim; 

% Preeipitaion Diameter vector (divider<D<mid) 
vector_size_2=diam_lim:resolution2:diam_mid; 

% Preeipitaion Diameter vector (mid<D<highest) 
vector_size_3=diam_mid:resolution3:diam_high; 

% Preeipitaion Diameter Vector 

vector_size=[vector_size_l,vector_size_2,vector_size_3]; 


% - 

% Create RCS vector for all Rain drops of Precipitation Vector. 
rcs=((kw^2)*(pU5)/(wavelength^4))*((vector_size*le-3).^6); 
%- 


% - 

% The Marshall Palmer approximation is used to create Preeipitaion Diameter 
% Distribution Vector (Drop Size Distribution N(D) Vector) 

% The number of drops are scaled down by a factor (Norm) to improve the 
% speed of the simulator 

% - 

lambda=4.1 * (rain_rate^(-0.21)); % [mm-1 ] 

number=8000*exp(-lambda*vector_size)/Norm; 

vectorres 1 =ones( 1 ,length(vector_size_l ))*resolutionl; 
vector_res2=ones(l,length(vector_size_2))*resolution2; 
vector_res3 =ones( 1, length(vector_size_3)) *resolution3; 

% delta(Diameter) dD 

vector_res=[vector_resl,vector_res2,vector_res3]; 

% the approximate "true" number of drops N(D)dD vector 
number=number.*vector_res; 

% - 


% - 

% create matrix T 
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T=[vector_size;rcs;number]; 


% - 

% Change matrix 'T' to matrix 'matrix T Round' 

% - 

% This part computes the total number of drops by utilizing the resolution 
% cell and the matrix T Round 

% The only difference between the matrix T Round and T, is the fact that 
% matrix T Round has the total number of drops on the third row and they 
% are rounded of to nearest integer value. 

% - 


matrix_T_Round=T; 

matrix_T_Round(3,:)=round(matrix_T_Round(3,:)*volume_box); 

% Keep track of errors due to round off 
diff=T (3, :)* volume_box-matrix_T_Round(3, 

error_rcs=sum(diff *matrix_T_Round(2,:))*Norm/volume_box*volume_resolution; 
error_Z=error_rcs/volume_resolution*wavelength^4/pi^5/abs(kw)*2/le-18; 


% - 

% Part III, Initial positioning of the drops within the volume 

% - 

% This part generates the initial position of the drops in the volume of 
% the box. 

% The positions are recorded in the fde "sizexx.maf taking into account 
% the azimuth angle 

% - 

% INPUT 
% - box 

% - matrix T Round 
% - elevation 
% - azimuth 

% - 

% OUTPUT 

% All fdes including the sizes. Each contains the positions of all drops 
% for a given diameter. The data recoreded as follows: 

% - radial distance r 
% - elevation phi 
% - azimuth theta 

% - angle with respect to beam center. This parameter will be 
% used to identify drops which are in the beam. 

% - 


% Randomly put the calculated number of drops of each diameter in the box 
number_size=length(matrix_T_Round( 1,:; 

for i=l :number_size 

% select number of drops of certain Diameter 
drops=matrix_T_Round(3 ,i); 

% create as many positions on the x axis as # drops per size 
elements_x=rand( 1,drops); 

% put them in the range of interest 
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% Random placement x for every drop 

position_x=box( l,l)+(box(l,3)-box(l,l))* elementsx; 

% create as many positionson the y axis as # drops per size 
elements_y=rand( 1,drops); 

% put them in the range of interest 
% Random placement y for every drop 

position_y=box(2,l)+(box(2,2)-box(2,l))*elements_y; 

% create as many positions on the z axis as # drops per size 
elements_z=rand( 1,drops); 

% put them in the range of interest 
% Random placement z for every drop 

position_z=box(3,1 )+(box(3,5)-box(3,1 ))*elements_z; 

% matrix of results 

% Creating a drop Matrix using spherical coordinates 

[theta,phi,r]=cart2sph(position_x,position_y,position_z); 

% Calculating the range, phi, theta and angle to every drop with 
% reference to phi O and theta O. Used to estimate weighted power return. 
vector_azimuth=ones( 1 ,drops) * ((azimuth *pi)/180); 
vector_elevation=ones(l,drops)*((elevation*pi)/180); 
[xf,yf,zf|=sph2cart(vector_azimuth,vector_elevation,r); 
range=sqrt((position_x-xf). ^2 .+(position_y-y f) .^2. +(position_z-zf) .^2); 

% approximation 
angle=atan(range. /r); 

% Saving the drop positions in size%d.mat 
% matrixpos 

matrixpos=[r;phi;theta;angle]; 
save(sprintf('size%d.mat',i),'matrixpos'); 
end 

% - 

% - 

% Part IV, Calculate Coherent Electric field return. 

% - 

% For every pulse the coherent Electric field return is calculatet. Between 
% pulses all drops are moved. 

% - 


% Run the simulation for specified number of pulses 

% Zero all vectors and matrices that are going to be used 

phase_vector=[]; 

sum_E=[]; 

weighted_E=[]; 

for p=l :number_pulses 

disp(sprintf('calculate the power-return of pulse %d\n',p)); 
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% compute the (weighted) I and Q return 

% - 

% computes the total number of drops 

% - 

% INPUT 
% - Matrix T 
% - wavelength 
% - angle_3dB 
% - distance 
% - coord vol res 
% - number_pulses 

% - 

% OUTPUT 

% - total E: the received complex Electric field) 

% - power: the power 
% - voltage: the voltage 

% - drops beam: the number of drops in the beam 

% - 

% Only the drops which are in the resolution cell are selected for the 
% computation 

% - 


number_size=length(matrix_T_Round( 1,:)); 

total_E=0; 

total_phase=0; 

new_E=0; 

for i=l:number_size 

% opening the file of registered data matrixpos. 
load(sprintf('size%d.mat',i)); 

% removing drops that are outside the resolution cell 

number_drops_per_size=length(matrixpos(l,:)); 

dropmatrix=[]; 

% Filling the drop matrix with all drops inside the beam 

% resolution volume 

for d=l :number_drops_per_size 

if ((matrixpos(4,d)<=(angle_3db*pi/180))&(matrixpos(l,d)>^-- 

coord_vol_res(l, 1 ))&(matrixpos( 1 ,d)<=coord_vol_res( 1,2))); 
dropmatrix=[dropmatrix,matrixpos(: ,d)]; 
end 
end 


% - 

% Convert to cartesian coordinates 
if isempty(dropmatrix)==0 

[position_x,position_y,position_z]=sph2cart(dropmatrix(3,:),... 
dropmatrix(2, :),dropmatrix( 1,:)); 

% Create comparison vectors 

vector_azimuth=ones(l ,length(dropmatrix(l, :)))*((azimuth*pi)/l 80); 
vector_elevation=ones( 1 ,length(dropmatrix( 1,:))) * ((ele vation*pi)/180); 
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% Create a phase veetor for the drops of eurrent Diameter 
r_phase=position_x. *u+position_y. *v+position_z. *w; 
phase=exp((2 *j *k). *r_phase); 

% Sum the eontributions to the E-field 

total_phase=total_phase+sum(phase); 

total_E=total_E+sum(sqrt(matrix_T_Round(2,i)).*phase); 


% - 

% eeate weighted power return, eoneidering range (r), theta, phi 
% with referenee to boresight and radar parameters 
E_weighted=sqrt(P_t)*G*wavelength/(sqrt(4*pi)^3)./dropmatrix(l,:)... 

. ^2. *(exp(-4 *log(2). * ((dropmatrix(3,: )-veetor_azimuth) .^2/... 
(2*angle_3db*pi/180)^2+((dropmatrix(2,:)-veetor_elevation)... 

. ^2/(2 *angle_3 db*pi/180)*2)))); 

new_E=new_E+sum(E_weighted. * (sqrt(matrix_T_Round(2,i)). *phase)); 

end 

end 

% - 


% - 

% This part simulates the movement of the drops and gives new positions. 

% - 

% INPUT 

% - matrix T Round 
% - PRT 
% - azimuth 
% - elevation 
% - rain_anglel 
% - rain_angle2 
% - diam lim 
% - resolution! 

% - resolutioiiZ 

% - wind veetor 

% - 

% OUTPUTS 

% Output files eontaining sizes with new positions 
% struetured the same way as in previous 

% there is also an output fde for speed (one for eaeh size) for the pulse 
% pair. 

% - 


number_size=length(matrix_T_Round( 1,:)); 

% adding a Gaussian random term to make the wind speed spread 
if wind_veetor( 1 )~=0 

new_wind_veetor( 1 )=wind_veetor( 1 )*normmd( 1 ,abs(spread/max(wind_veetor( 1)))); 
else 

new_wind_veetor( 1 )=0; 
end 

if wind_veetor(2)~=0 

new_wind_veetor(2)=wind_veetor(2)*normmd(l,abs(spread/max(wind_veetor(2)))); 

else 

new_wind_veetor(2)=0; 
end 
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if wind_vector(3)~=0 

new_wind_vector(3)=wind_vector(3)*normmd(l,abs(spread/max(wind_vector(3)))); 

else 

new_wind_veetor(3)=0; 

end 

for i=l:number_size 

% opening of data fde registered in matrixpos. 
load(sprintf('size%d.mat',i)); 

% total number of drops 
drops=matrix_T_Round(3 ,i); 

% transforming to reetangular eoordinates 
[x,y,z]=sph2eart(matrixpos(3,:),matrixpos(2,:),matrixpos(l,:)); 

diameter=matrix_T_Round( 1 ,i); 


% - 

% Drop fall speed depends on drop diameter in the z axis, 

% Atlas-Ulbrieh approximation is used 

% Approximation is valid in the diameter range 5*e-4, 5*e-3 
%- 


wtmax=386.6*((matrix_T_Round(l,i)*0.001)^0.67); 

veloeityz=(wtmax*(ones(l,drops)))*eos((rain_anglel*pi)/180); 

veloeityx=(wtmax*(ones(l,drops)))*sin((rain_anglel*pi)/180)*... 

eos((rain_angle2*pi)/180); 

veloeityy=(wtmax*(ones(l,drops)))*sin((rain_anglel*pi)/180)*... 
sin((rain_angle2*pi)/l 80); 


newest_wind_veetor=[]; 


% adding Gaussian random term to make the wind speed 
% spread between drops 
if new_wind_veetor(l)~=0 

newest_wind_veetor( 1, :)=new_wind_veetor( 1) *... 

normmd(l,abs(spread_2/max(new_wind_veetor(l))),l,drops); 

else 

newest_wind_veetor( 1,: )=0; 
end 

if new_wind_veetor(2)~=0 

newest_wind_veetor(2,:)=new_wind_veetor(2)*... 

normmd(l,abs(spread_2/max(new_wind_veetor(2))),l,drops); 

else 

newest_wind_veetor(2,: )=0; 
end 

if new_wind_veetor(3)~=0 

newest_wind_veetor(3, :)=new_wind_veetor(3) *... 

normmd(l,abs(spread_2/max(new_wind_veetor(3))),l,drops); 

else 

newest_wind_veetor(3,: )=0; 
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end 


% If spread within pulse 
veloeityxt=newest_wind_veetor(l,:)-veloeityx; 
veloeityyt=newest_wind_veetor(2,:)-veloeityy; 
veloeityzt=newest_wind_veetor(3,:)-veloeityz; 


% speed registration 

m_veloeity=[veloeityxt;veloeityyt;veloeityzt]; 

save(sprintf('veloeity%d.mat',i),'m_velocity'); 

deplaeementx=veloeityxt*PRT; 
deplaeementy=velocityyt*PRT; 
deplaeementz=veloeityzt*PRT; 

position_x=x+deplaeementx; 

position_y=y+deplacementy; 

position_z=z+deplaeementz; 

% transformation into spherieal eoordinates 

[theta,phi,r]=eart2sph(position_x,position_y,position_z); 

% eomputation of new angles 
veetor_azimuth=ones( 1 ,drops) *(azimuth*pi)/180; 
veetor_elevation=ones(l,drops)*(elevation*pi)/180; 
[xf,yf,zf|=sph2eart(veetor_azimuth,veetor_elevation,r); 

range=sqrt((position_x-xf).^2.+(position_y-yf).^2.+(position_z-zf).^2); 
% approximation 
angle=atan(range ./r); 

% matrixpos 

matrixpos=[r;phi;theta;angle]; 

% save to file, matrixpos 

s=sprintf('size%d. mat',i); 
save(s, 'matrixpos'); 


end 

%■ 


% add to E-field and phase veetor 
sum_E=[sum_E,total_E]; 
weighted_E=[weighted_E,new_E]; 
phase_veetor=[phase_veetor,ANGLE(total_phase)]; 

end 


%■ 


% - 

% Power return for every pulse ineluding error 

running_Z=(Norm*abs(sum_E).^2)*wavelength^4/pi^5/abs(kw)^2/volume_resolution... 

/le-18+error_Z; 
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power_retum=Norm*abs(weighted_E). ^2; 


% [W] 


% calculate average Power and Z 
Z_avg=sum(running_Z)/p; 

Z_avg_plot=ones( 1 ,p) *Z_avg; 

power_retum_avg=sum(Norm*abs(weighted_E).^2)/p; % [W] 

power_retum_avg_plot=ones(l,p)*power_retum_avg; 

% calculate doppler frequency 
phase_vector_u=unwrap(phase_vector); 

doppler_fq=([phase_vector_u(2 :p)-phase_vector_u( 1 :p-1 )])/2/pi/PRT; 
number_vector=l :p; 


% - 

% Calculate references 
% Integral form of Z 
11 = 0 : 100 / 1000 : 100 ; 

zz=11.^6.*8000.*exp(-4.1.*raln_rate *(-.21).*ll); 
real_Z=ones( 1 ,number_pulses)* 10*log 10(trapz(ll,zz)); 

% reflectivity_real=ones(l,number_pulses)*10*logl0(pC5/wavelength^4*... 

% volume_resolutlon*abs(kw)^2*real_Z* le-18); 

% Adding the error to the welghted E estimations 

error_power=error_rcs*P_t*G'^2*wavelength*2/pC2*(2*angle_3db*pl/l 80)^2*... 
(c*pulsewldth*le-6)/1024/log(2)/pC2/(dlstance*le3)^2/volume_resolutlon; 

% Evaluate average Power return to make Z estimation Including error 
Z_estlmatlon=power_retum_avg/P_t/G^2/wavelength^2*... 

pC2/(2*angle_3db*pl/180)^2/(c*pulsewldth*le-6)*1024*log(2)*(dlstance*le3)^2/... 
((kw''2)*(pC5)/(wavelength^4))/le-18; 

Z_estlmatlon_plot=ones( 1 ,number_pulses)* 10*log 10(Z_estlmatlon+error_Z); 

% calculate the EFT of the weighted electric field return (1024 points) to 
% get the power spectal estimate 
NN=1024; % number of FFT points 

spectral_estlmate=abs(fft(welghted_E-mean(welghted_E),NN)).^2; 

spectral_estlmate_shlfted=fftshlft(spectral_estlmate); 

% - 

% Part V, Plot of results 

% - 


figure(l); 

subplot(2,l,l); 

plot(number_vector(l:p-l),doppler_fq(l:p-l),'-x',number_vector(l:p),... 
ones(l,p)*mean(doppler_fq),'—r'); 

tltle(['Doppler frequency for all pulse pairs for PRF= ',num2sfr(l/PRT),' Hz']); 

xlabel('Pulse pair'); 

ylabel('f_d'); 

legend('Doppler frequency per pulse pair','Average Doppler estimation') 
grid on 


subplot(2,l,2); 

v_r=doppler_fq* wavelength/2; 
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hist(v_r,15) 
hold on 


[ff,xx] = ksdensity(v_r); 
plot(xx,ff/max(ff)*max(hist(v_r,15))) 
tide('Radial Velocity Histogram with Gaussian Fit'); 
xlabel('Radial velocity [m/s]'); 
ylabel('Frequency'); 
grid on 

figure (2) 
subplot(2,l,l) 

% also converting to dBm 

plot(number_vector,10*logl0(power_retum+error_power)+30,'-x',... 

number_vector,10*logl0(power_retum_avg_plot+error_power)+30,'—r'); 
title('Power Return for all pulses'); 
xlabel('Sample time'); 
ylabel('Power [dBm]'); 

legend('Power return', 'Average Power return',0) 
grid on 

subplot(2,l,2); 

plot(number_vector*1000*PRT,10*logl0(abs(mnning_Z)),'-xb',number_vector... 
*1000*PRT,10*logl0(abs(Z_avg_plot)),'-.r',number_vector*1000*... 
PRT,real_Z,'-k',number_vector*1000*PRT,Z_estimation_plot,':b'); 
title('Estimated Z per pulse'); 
xlabel('Sample time'); 
ylabel('Z [dBZ]'); 

legend('Simulation Z, actual', 'Average Z','Z=\int D^6N(D)dD',... 

'Estimated Z, Zeroth moment',0) 
hbne = findobj(gca,'Type','line'); 
set(hline,'LineWidth',2) 
grid on 


figure (3) 
subplot(2,l,l) 

f=abs(xcorr((weighted_E),'coeff)); 
tau_ts=ones( 1 ,length(f))*f((lengtli(f)+1 )/2-1); 

plot(l:length(f),abs(xcorr((weighted_E),'coeff)),'b',l:length(f),tau_ts,'-r'); 
title('Retum Voltage Correlation'); 
xlabel('Number of pulses'); 
ylabel('Normalized correlation value'); 

legend(['Correlation function'],['\rho(\tau=T_s)=',num2str(tau_ts( 1))]) 
grid on 

subplot(2,l,2) 

% Prepare for ticks so every 100 Hz is marked 
n=0:l/PRT/100; 

N_i=100*(NN-l)*n*PRT; 
alfa=zeros(l,l/PRT/100+l); 
vic=(length(alfa)-1 )/2; 
for beta=l :length(alfa) 
alfa( 1 ,beta)=-vic* 100; 
vic=vic-l; 
end 
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avg_fq=mean(doppler_fq) 
std_fq=std(doppler_fq) 
gauss_vector=min(alfa): 10: max(alfa); 
uu=length(gauss_vector); 

[value l,mdexl]=max(abs(spectral_estimate_shifted)); 
avg_fq_fft=(index 1 -NN/2)/PRT/NN; 

plot_gauss_vector=linspace(0,NN,uu); 

plot_nom_fq=normpdf(gauss_vector,avg_fq,std_fq/2); 

plot_nom_fft=nompdf(gauss_vector,avg_fq_fft,std_fq/2); 

plot(spectral_estimate_shifted/max(spectral_estimate_shifted)) 
hold on 

plot(fftshift(abs(fft(xcorr((weighted_E),'coeff),1024))/... 

max(abs(fft(xcorr((weighted_E),'coeff),1024)))),'-k') 
title('Power Spectral Estimate'); 
xlabel('Doppler Frequancy [Flz]'); 
ylabel('S|fl, normalized'); 
xlim([0 NN]); 
set(gca,'xtick', N_i) 
set(gca,'xticklaber, num2str(alfa')); 
grid on 

2) WR_dual_prf 

% Weather Radar Signal Simulator, developed version 
% Two PRFs enabling pulse pair comparison 
% Ulf Schroder 
% August, 2005 

% - 


% Main program for The NPS/Marielle Gosset Weather Radar Signal Simulator. 


% - 

% DESCRIPTION 
% Part I 

% From Radar Data: 

% - Generate Beam Resolution cell approximation 

% - Generate Cube adding a margine to the Beam Resolution cell to allow for 
% the rain drops to fall trough the through beam over time. 

% 

% Part II 

% From Rain Parameters: 

% Generating the fundamental T Matrix concisting of 
% - Diameter Vector and Number of drops in each diameter interval 
% using the Marshall-Palmer exponential approximation. 

% - the backscatter RCSs for each Diameter size 
% - the total number of drops for each size per m''3 

% Part III 

% Randomly place all drops in cube 
% Part IV 

% For every pulse the coherent Electric field return is calculated. Between 
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% pulses all drops are moved. 


% Part V 

% Caleulate and plot results 

% - 

elear all 
elose all 
ele 

% Input Data 

e=3e8; % Speed of light [m/s] 

% Radar Data 

P_t= 10e3; % Transmitted power 

G_dB=20; % Maximum gain in dB 

G=10^(G_dB/10); % Maximum gain 

angle_3db=0.20; % 3 dB Beam angle [degrees], (half half power beamwidth) 

pulsewidth=0.20; % Pulsewidth [mieroseeond] 

elevation=3; % Elevation angle [degrees] 

margine=5; % margin for resolution eell[m] 

azimuth=0; % Azimuth angle [rad] 

number_pulses=40; % Number of pulses 

PRT_l=l/2000; % Doppler Pulse Repetition Time [seeonds] 

PRT_2=l/200; % Avg Power Pulse Repetition Time [seeonds] 

distanee=0.5; % Range from the radar to the resolution eell [Km] 
wavelength=0.1; % Wavelength [m] 

k=2*pi/wavelength; 

% Rain parameters 

diam_low=0.05; % Lowest diameter of Preeipitation [mm] 

diam_hm=l; % Divider between low and mid [mm] 

diam_mid=5; % Divider between mid and high [mm] 

diam_high=15; % Highest diameter of Preeipitation [mm] 

resolutionl=0.08; % Lower diameter resolution [mm] 
resolution2=0.3; % Higher diameter resolution [mm] 
resolution3=0.8; % Highest diameter resolution [mm] 

sort_approx=l; % RCS approximation model; [l]=Rayleigh 

kw=0.93; % abs(K)*2 approx for water 

rain_rate=50; % Rain Rate [mm/h] 
rain angleDO; % degrees Rain Angle (Phi) 
rain_angle2=30; % degrees Rain Angle (Theta) 
wind_veetor=[-10;0;0]; % Wind speed veetor (x,y,z) [m/s] 
spread=l; % Wind speed spread [m/s], (all direetions) 

spread_2=4; % Drop veloeity spread [m/s], 

Norm=100; % Normalization eoeffieient 

% Direetional eosines 

u=sin(pi/2-elevation*pi/180)*eos(azimuth*pi/180); 
v=sin(pi/2-elevation*pi/180) * sin(azimuth*pi/l 80); 
w=eos(pi/2-elevation*pi/180); 

% - 

% Part I, Generate the resolution volume 

% - 

% Approximation of the resolution eell of the pulse eireular box, 
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% for any azimuth angle 

% - 

% INPUT 
% - distance 
% - angle_3dB 
% - pulsewidth 
% - elevation 
% - margine 
% - azimuth 

% - 

% OUTPUT 

% - volume resolution: computed volume of the resolution cell 
% - volume box: volume of the box where the drops are put 
% - box: the 8 rectangular coordinates of the box 
% - coord vol res: the coordinates of the resolution cell where the drops 
% are taken. 

% - 


% - 

% Creating the circular beam resolution cell 
range_res=(c*pulsewidth* le-6)/2; 

% Range resolution (ct/2) [m] 

% small angle approximation allows for following 
angle_rad=(angle_3 db *pi)/180; 

% 3 dB beamwidth [rad] 
radius 1=(tan(angle_rad)* (distance *1000)); 

% Radius of beam at R=distance 
radius2=(tan(angle_rad)*((distance* 1000)+range_res)); ... 

% Radius of beam at R=distance + range resolution 

surface 1 =pi *radius U2; 
surface2=pi *radius2^2; 

% volume of cone = base surface times height for all three 
volumecone 1 =(surface 1 *(distance* 1000))/3; 
volume_cone2=(surface2*((distance*1000)+range_res))/3; 
volume_resolution=(volume_cone2-volume_conel); % Beam Cone Volume 

% computation of the resolution cell coordinates, spherical coordinates 
% elevation in radians 

elevation_rad=(elevation*pi)/180; % Elevation [rad] 

azimuth_rad=(azimuth*pi)/180; % Azimuth [rad] 

% Matrix used when comparing wheather a drop is inside or outside 
% beamwidth 

coord_vol_res=[distance* 1000,(distance* 1000)+range_res;elevation_rad-angle_rad,... 

elevation_rad+angle_rad;azimuth_rad-angle_rad,azimuth_rad+angle_rad]; 

% - 


% - 

% Creating the box 

% computation of the box and its coordinates 

coinl=[azimuth_rad+angle_rad,elevation_rad+angle_rad,distance* 1000]; 
coin2=[azimuth_rad-angle_rad,elevation_rad+angle_rad,distance* 1000]; 
coin3=[azimuth_rad+angle_rad,elevation_rad-angle_rad,distance *1000]; 
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coin4=[azimuth_rad-angle_rad,elevation_rad-angle_rad,distance* 1000]; 
coin5=[azimuth_rad+angle_rad,elevation_rad+angle_rad,distance* 1000+range_res]; 
coin6=[azimuth_rad-angle_rad,elevation_rad+angle_rad,distance*1000+range_res]; 
coin7=[azimuth_rad+angle_rad,elevation_rad-angle_rad,distance*1000+range_res]; 
coin8=[azimuth_rad-angle_rad,elevation_rad-angle_rad,distance* 1000+range_res]; 

[xl,yl,zl]=sph2cart(coinl(l),coinl(2),coinl(3)); 

[x2,y2,z2]=sph2cart(coin2(l),coin2(2),coin2(3)); 

[x3 ,y3 ,z3]=sph2cart(coin3 (l),coin3 (2),coin3 (3)); 
[x4,y4,z4]=sph2cart(coin4(l),coin4(2),coin4(3)); 
[x5,y5,z5]=sph2cart(coin5(l),coin5(2),coin5(3)); 
[x6,y6,z6]=sph2cart(coin6(l),coin6(2),coin6(3)); 

[x7,y7,z7]=sph2cart(coin7 (l),coin7 (2),coin7 (3)); 
[x8,y8,z8]=sph2cart(coin8(l),coin8(2),coin8(3)); 

xmax=max([xl,x2,x3,x4,x5,x6,x7,x8]); 

xmin=min([xl,x2,x3,x4,x5,x6,x7,x8]); 

ymax=max([yl,y2,y3,y4,y5,y6,y7,y8]); 

ymin=min([yl,y2,y3,y4,y5,y6,y7,y8]); 

zmax=max([zl,z2,z3,z4,z5,z6,z7,z8]); 

zmin=min([zl,z2,z3,z4,z5,z6,z7,z8]); 

% the box with its margin 

coin_l=[xmin-margine,ymin-margine,zmin-margine]; 

coin_2=[xmin-margine,ymax+margine,zmin-margine]; 

coin_3=[xmax+margine,ymax+margine,zmin-margine]; 

coin_4=[xmax+margine,ymin-margine,zmin-margine]; 

coin_5=[xmin-margine,ymin-margine,zmax+margine]; 

coin_6=[xmin-margine,ymax+margine,zmax+margine]; 

coin_7=[xmax+margine,ymax+margine,zmax+margine]; 

coin_8=[xmax+margine,ymin-margine,zmax+margine]; 


box=[coin r,coin_2',coin_3',coin_4',coin_5',coin_6',coin_7',coin_8']; 

% calculate the volume of the box 

volume_box=abs((xmax+margine-(xmin-margine))*(ymax+margine-(ymin-margine))*(... 

zmax+margine-(zmin-margine))); 


cube=[coinr,coin_2',coin_3',coin_4',coinr,coin_5',coin_6',coin_2',... 

coin_6',coin_7',coin_3',coin_7',coin_8',coin_4',coin_8',coin_5']; 

% - 


% - 

% Part 11, Generate the fundamental T Matrix 

% - 

% This part produces a matrix of all drops that contribute to the total 
% reflectivity of the resolution cell 

% - 

% First all sizes are created, with two vectors of different size resolution. 

% The size vectors represents all diameters. Next the backscatter RCS for each 
% size is computed. 

% Using the Marshall-Palmer exponential approximation, the number of drops of 
% each diameter is computed. 

% - 
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% INPUT 
% - diam low 
% - diam lim 
% - diam high 
% - resolution! 

% - resolution2 
% - sort approx 
% - wavelength 
% - Kw 
% - rain_rate 

% - 

% OUTPUT 
% - number drops 
% Matrix T: 

% - first row all sizes 
% - second row all the backscatter RCSs 
% - third row the total number of drops for each size per m^3 
% the columns start with the smallest size 

% - 

% Create dropsize vector 

% Precipitaion Diameter vector (lowest<D<divider) 
vector_size_ 1 =diam_low resolution 1 :diam_lim; 

% Precipitaion Diameter vector (divider<D<mid) 
vector_size_2=diam_lim:resolution2:diam_mid; 

% Precipitaion Diameter vector (mid<D<highest) 
vector_size_3=diam_mid:resolution3:diam_high; 

% Precipitaion Diameter Vector 

vector_size=[vector_size_l,vector_size_2,vector_size_3]; 


% - 

% Create RCS vector for all Rain drops of Precipitation Vector. 
rcs=((kw''2)*(pU5)/(wavelength^4))*((vector_size*le-3).^6); 
%- 


% - 

% The Marshall Palmer approximation is used to create Precipitaion Diameter 
% Distribution Vector (Drop Size Distribution N(D) Vector) 

% The number of drops are scaled down by a factor (Norm) to improve the 
% speed of the simulator 

% - 

lambda=4.1 * (rain_rate^(-0.21)); % [mm-1 ] 

number=8000*exp(-lambda*vector_size)/Norm; 

vectorres 1 =ones( 1 ,length(vector_size_l ))*resolutionl; 
vector_res2=ones(l,length(vector_size_2))*resolutiorL2; 
vector_res3=ones(l,length(vector_size_3))*resolution3; 

% delta(Diameter) dD 

vector_res=[vector_resl,vector_res2,vector_res3]; 

% the approximate "true" number of drops N(D)dD vector 
number=number.*vector_res; 
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%■ 


% - 

% create matrix T 
T=[vector_size;rcs;number]; 


% - 

% Change matrix 'T' to matrix 'matrix T Round' 

% - 

% This part computes the total number of drops by utilizing the resolution 
% cell and the matrix T Round 

% The only difference between the matrix T Round and T, is the fact that 
% matrix T Round has the total number of drops on the third row and they 
% are rounded off to nearest integer value. 

% - 


matrix_T_Round=T; 

matrix_T_Round(3,:)=round(matrix_T_Round(3,:)*volume_box); 

% Keep track of errors due to round off 
diff=T (3, :)* volume_box-matrix_T_Round(3, 

error_rcs=sum(diff *matrix_T_Round(2,:))*Norm/volume_box*volume_resolution; 
error_Z=error_rcs/volume_resolution*wavelength^4/pi^5/abs(kw)*2/le-18; 


% - 

% Part III, Initial positioning of the drops within the volume 

% - 

% This part generates the initial position of the drops in the volume of 
% the box. 

% The positions are recorded in the fde "sizexx.maf taking into account 
% the azimuth angle 

% - 

% INPUT 
% - box 

% - matrix T Round 
% - elevation 
% - azimuth 

% - 

% OUTPUT 

% All files including the sizes. Each contains the positions of all drops 
% for a given diameter. The data recoreded as follows: 

% - radial distance r 
% - elevation phi 
% - azimuth theta 

% - angle with respect to beam center. This parameter will be 
% used to identify drops which are in the beam. 

% - 


% Randomly put the calculated number of drops of each diameter in the box 
number_size=length(matrix_T_Round( 1,:; 

for i=l inumber size 

% select number of drops of certain Diameter 
drops=matrix_T_Round(3 ,i); 
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% create as many positions on the x axis as # drops per size 
elements_x=rand( 1,drops); 

% put them in the range of interest 
% Random placement x for every drop 

position_x=box( 1,1 )+(box( 1,3)-box(l, 1 ))*elements_x; 

% create as many positionson the y axis as # drops per size 
elements_y=rand( 1,drops); 

% put them in the range of interest 
% Random placement y for every drop 

position_y=box(2,l)+(box(2,2)-box(2,l))*elements_y; 

% create as many positions on the z axis as # drops per size 
elements_z=rand( 1,drops); 

% put them in the range of interest 
% Random placement z for every drop 

position_z=box(3,1 )+(box(3,5)-box(3,1 ))*elements_z; 

% matrix of results 

% Creating a drop Matrix using spherical coordinates 

[theta,phi,r]=cart2sph(position_x,position_y,position_z); 

% Calculating the range, phi, theta and angle to every drop with 
% reference to phi O and theta O. Used to estimate weighted power return. 
vector_azimuth=ones( 1 ,drops) * ((azimuth *pi)/180); 
vector_elevation=ones(l,drops)*((elevation*pi)/180); 
[xf,yf,zf|=sph2cart(vector_azimuth,vector_elevation,r); 
range=sqrt((position_x-xf) .^2 .+(position_y-y f) .^2 .+(position_z-zf) .^2); 

% approximation 
angle=atan(range. /r); 

% Saving the drop positions in size%d.mat 
% matrixpos 

matrixpos=[r;phi;theta;angle]; 
save(sprintf('size%d.mat',i),'matrixpos'); 
end 

% - 


% - 

% Part IV, Calculate Coherent Electric field return. 
%- 


% For every pulse the coherent Electric field return is calculatet. Between 
% pulses all drops are moved. 


%■ 


% Run the simulation for specified number of pulses 

% Zero all vectors and matrices that are going to be used 

powerT=0; 

vector_power=[]; 
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phase_vector=zeros(l,2*number_pulses); 

imaginary_retum=zeros(l,2*number_pulses); 

sum_E=[]; 

weigbted_E=[]; 

for p=l :number_pulses 

disp(spriiitf('calculate tbe power-rebim of pulse %d\n',p)); 

% compute tbe (weighted) I and Q return 

% - 

% computes tbe total number of drops 

% - 

% INPUT 
% - Matrix T 
% - wavelength 
% - angle_3dB 
% - distance 
% - coord vol res 
% - number_pulses 

% - 

% OUTPUT 

% - total E: the received complex Electric field) 

% - power: the power 
% - voltage: the voltage 

% - drops beam: the number of drops in the beam 

% - 

% Only the drops which are in the resolution cell are selected for the 
% computation 

% - 


number_size=length(matrix_T_Round( 1,:)); 

total_E=0; 

total_phase=0; 

new_E=0; 

for i=l :number_size 

% opening the file of registered data matrixpos. 
load(sprintf('size%d.mat'd)); 

% removing drops that are outside the resolution cell 

number_drops_per_size=length(matrixpos(l,:)); 

dropmatrix=[]; 

% Filling the drop matrix with all drops inside the beam 

% resolution volume 

for d=l :number_drops_per_size 

if ((matrixpos(4,d)<=(angle_3db*pi/180))&(matrixpos(l,d)>^-- 

coord_vol_res(l, 1 ))&(matrixpos( 1 ,d)<=coord_vol_res( 1,2))); 
dropmatrix=[dropmatrix,matrixpos(:,d)]; 
end 
end 


% - 

% Convert to cartesian coordinates 
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if isempty(dropmatrix)==0 

[position_x,position_y,positioii_z]=sph2cart(dropmatrix(3, 
dropmatrix(2, :),dropmatrix( 1, 

% Create comparison vectors 

vector_azimuth=ones(l ,length(dropmatrix(l, :)))*((azimuth*pi)/l 80); 
vector_elevation=ones( 1 ,length(dropmatrix( 1,:))) * ((ele vation*pi)/180); 

% Create a phase vector for the drops of current Diameter 
r_phase=position_x.*u+position_y.*v+position_z.*w; 
phase=exp((2 *j *k). *r_phase); 

% Sum the contributions to the E-field 

total_phase=total_phase+sum(phase); 

total_E=total_E+sum(sqrt(matrix_T_Round(2,i)).*phase); 


% - 

% ceate weighted power return, concidering range (r), theta, phi 
% with reference to boresight and radar parameters 

E_weighted=sqrt(P_t)*G*wavelength/(sqrt(4*pi)^3)./dropmatrix(l,:).''2.*... 
(exp(-4*log(2).*((dropmatrix(3,:)-vector_azimuth).*2/... 
(2*angle_3db*pi/180)^2+((dropmatrix(2,:)-vector_elevation) *2/... 
(2*angle_3db*pi/180)''2)))); 

new_E=new_E+sum(E_weighted. * (sqrt(matrix_T_Round(2,i)). *phase)); 

end 

end 

phase_vector(2*p-l)=ANGLE(total_phase); 
imaginary_retum(2 *p-1 )=total_phase; 


%■ 


% - 

% This part simulates the movement of the drops and gives new positions. 
% First using the Doppler PRT 

% - 

% INPUT 

% - matrix T Round 
% - PRT 
% - azimuth 
% - elevation 
% - rain_anglel 
% - rain_angle2 
% - diam lim 
% - resolution! 

% - resolutiorrZ 

% - wind vector 

% - 

% OUTPUTS 

% Output files containing sizes with new positions 
% structured the same way as in previous 

% there is also an output fde for speed (one for each size) for the pulse 
% pair. 

% - 
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number_size=length(matrix_T_Round( 1,:; 

% adding a Gaussian random term to make the wind speed spread 
if wind_vector( 1 )~=0 

new_wind_vector( 1 )=wind_vector( 1 )*normmd( 1 ,abs(spread/max(wind_veetor( 1)))); 
else 

new_wind_vector( 1 )=0; 
end 

if wind_vector(2)~=0 

new_wind_vector(2)=wind_vector(2)*normmd(l,abs(spread/max(wind_vector(2)))); 

else 

new_wind_vector(2)=0; 
end 

if wind_vector(3)~=0 

new_wind_vector(3)=wind_vector(3)*normmd(l,abs(spread/max(wind_vector(3)))); 

else 

new_wind_vector(3)=0; 

end 


for i=l mumber size 

% opening of data fde registered in matrixpos. 
load(sprintf('size%d.mat',i)); 

% total number of drops 
drops=matrix_T_Round(3 ,i); 

% transforming to reetangular coordinates 
[x,y,z]=sph2cart(matrixpos(3,:),matrixpos(2,:),matrixpos(l,:)); 

diameter=matrix_T_Round( 1 ,i); 


% - 

% Drop fall speed depends on drop diameter in the z axis, 

% Atlas-Ulbrich approximation is used 

% Approximation is valid in the diameter range 5*e-4, 5*e-3 
%- 


wtmax=386.6*((matrix_T_Round(l,i)*0.001)^0.67); 

velocityz=(wtmax*(ones(l,drops)))*cos((rain_anglel*pi)/180); 

velocityx=(wtmax*(ones(l,drops)))*sin((rain_anglel*pi)/180)*... 

cos((rain_angle2*pi)/180); 

velocityy=(wtmax*(ones(l,drops)))*sin((rain_anglel*pi)/180)*... 
sin((rain_angle2*pi)/l 80); 

newest_wind_vector=[]; 


% adding Gaussian random term to make the wind speed 
% spread between different drops sizes 
if new_wind_vector(l)~=0 

newest_wind_vector( 1, :)=new_wind_vector( 1) *normmd( I,... 
abs(spread_2/max(new_wind_vector(l))),l,drops); 

else 
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newest_wind_vector( 1,: )=0; 
end 

if new_wind_vector(2)~=0 

newest_wind_vector(2,:)=new_wind_vector(2)*normmd(l,... 
abs(spread_2/max(new_wind_vector(2))),l,drops); 

else 

newest_wind_vector(2,: )=0; 
end 

if new_wind_veetor(3)~=0 

newest_wind_vector(3, :)=new_wind_vector(3) *normmd( I,... 
abs(spread_2/max(new_wind_vector(3))),l,drops); 

else 

newest_wind_veetor(3,: )=0; 
end 

velocityxt=newest_wind_vector(l,:)-velocityx; 
velocityyt=newest_wind_vector(2,:)-velocityy; 
velocity zt=newest_wind_vector(3,:)-velocityz; 

% speed registration 

m_velocity=[velocityxt;velocityyt;velocityzt]; 

save(sprintf('velocity%d.mat',i),'m_velocity'); 

deplacementx=velocityxt*PRT_l; 
deplacementy=velocityyt*PRT_l; 
deplacementz=velocityzt*PRT_l; 

position_x=x+deplacementx; 

position_y=y+deplacementy; 

position_z=z+deplacementz; 

% transformation into spherical coordinates 

[theta,phi,r]=cart2sph(position_x,position_y,position_z); 

% computation of new angles 
vector_azimuth=ones( 1 ,drops) *(azimuth*pi)/180; 
vector_elevation=ones(l,drops)*(elevation*pi)/180; 
[xf,yf,zf|=sph2cart(vector_azimuth,vector_elevation,r); 

range=sqrt((position_x-xf). ^2. +(position_y-y f) .^2. +(position_z-zf). ^2); 
% approximation 
angle=atan(range ./r); 

% matrixpos 

matrixpos=[r;phi;theta;angle]; 

% save to file, matrixpos 

s=sprintf('size%d. mat',i); 
save(s, 'matrixpos'); 


end 


% - 

% Calculate phase for pulse pair 
total_phase=0; 
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for i=l:number size 


% opening the file of registered data matrixpos. 
load(sprintf('size%d.mat',i)); 

% removing drops that are outside the resolution eell 

number_drops_per_size=length(matrixpos(l,:)); 

dropmatrix=[]; 

% Filling the drop matrix with all drops inside the beam 

% resolution volume 

for d=l :number_drops_per_size 

if ((matrixpos(4,d)<=(angle_3db*pi/180))&(matrixpos(l,d)>^-- 

eoord_vol_res(l, 1 ))&(matrixpos( 1 ,d)<=eoord_vol_res( 1,2))); 
dropmatrix=[dropmatrix,matrixpos(: ,d)]; 
end 
end 


% - 

% Convert to eartesian eoordinates 
if isempty(dropmatrix)==0 

[position_x,position_y,position_z]=sph2eart(dropmatrix(3, 

dropmatrix(2,:),dropmatrix(l,:)); 

% Create eomparison veetors 

veetor_azimuth=ones(l ,length(dropmatrix(l, :)))*((azimuth*pi)/l 80); 
veetor_elevation=ones( 1 ,length(dropmatrix( 1,:))) * ((elevation*pi)/180); 


% Create a phase veetor for the drops of eurrent Diameter 
r_phase=position_x.*u+position_y.*v+position_z.*w; 
phase=exp((2 *j *k). *r_phase); 

% Sum the eontributions to the E-field 
total_phase=total_phase+sum(phase); 
end 
end 

phase_veetor(2 *p)=ANGLE(total_phase); 
imaginary_retum(2*p)=total_phase; 


% - 

% This part simulates the movement of the drops and gives new positions 
% using the Avg Power PRT 

% - 


number_size=length(matrix_T_Round( 1,:)); 


% adding a Gaussian random term to make the wind speed spread 
if wind_veetor( 1 )~=0 

new_wind_veetor( 1 )=wind_veetor( 1 )*normmd( 1 ,abs(spread/max(wind_veetor( 1)))); 
else 

new_wind_veetor( 1 )=0; 
end 
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if wind_vector(2)~=0 

new_wind_vector(2)=wind_vector(2)*normmd(l,abs(spread/max(wind_vector(2)))); 

else 

new_wind_veetor(2)=0; 
end 

if wind_veetor(3)~=0 

new_wind_veetor(3)=wind_veetor(3)*normmd(l,abs(spread/max(wind_veetor(3)))); 

else 

new_wind_veetor(3)=0; 

end 


for i=l :number_size 

% opening of data fde registered in matrixpos. 
load(sprintf('size%d.mat',i)); 

% total number of drops 
drops=matrix_T_Round(3 ,i); 

% transforming to reetangular eoordinates 
[x,y,z]=sph2eart(matrixpos(3,:),matrixpos(2,:),matrixpos(l,:)); 

diameter=matrix_T_Round( 1 ,i); 


% - 

% Drop fall speed depends on drop diameter in the z axis, 

% Atlas-Ulbrieh approximation is used 

% Approximation is valid in the diameter range 5*e-4, 5*e-3 
%- 


wtmax=386.6*((matrix_T_Round(l,i)*0.001)^0.67); 

veloeityz=(wtmax*(ones(l,drops)))*eos((rain_anglel*pi)/180); 
veloeityx=(wtmax* (ones( 1,drops))) * sin((rain_anglel*pi)/180)*... 
eos((rain_angle2*pi)/180); 

veloeityy=(wtmax*(ones(l,drops)))*sin((rain_anglel*pi)/180)*... 
sin((rain_angle2*pi)/l 80); 

newest_wind_veetor=[]; 


% adding Gaussian random term to make the wind speed 
% spread between different drops sizes 
if new_wind_veetor(l)~=0 

newest_wind_veetor( 1, :)=new_wind_veetor( 1) *normmd( 1,... 
abs(spread_2/max(new_wind_veetor(l))),l,drops); 

else 

newest_wind_veetor( 1,: )=0; 
end 

if new_wind_veetor(2)~=0 

newest_wind_veetor(2,:)=new_wind_veetor(2)*normmd(l,... 
abs(spread_2/max(new_wind_veetor(2))),l,drops); 

else 

newest_wind_veetor(2,: )=0; 
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end 

if new_wind_vector(3)~=0 

newest_wind_vector(3, :)=new_wind_vector(3) *normmd( I,... 
abs(spread_2/max(new_wind_vector(3))),l,drops); 

else 

newest_wind_vector(3,: )=0; 
end 

velocityxt=newest_wind_vector(l,:)+velocityx; 
velocityyt=newest_wind_vector(2,:)+velocityy; 
velocity zt=newest_wind_vector(3,:)+velocityz; 

% speed registration 

m_velocity=[velocityxt;velocityyt;velocityzt]; 

save(sprintf('velocity%d.mat',i),'m_velocity'); 

deplacementx=velocityxt*(PRT_2-PRT_l); 
deplacementy=velocityyt*(PRT_2-PRT_l); 
deplacementz=velocityzt* (PRT_2-PRT_ 1); 

position_x=x+deplacementx; 

position_y=y+deplacementy; 

position_z=z+deplacementz; 

% transformation into spherical coordinates 

[theta,phi,r]=cart2sph(position_x,position_y,position_z); 

% computation of new angles 
vector_azimuth=ones( 1 ,drops) *(azimuth*pi)/180; 
vector_elevation=ones(l,drops)*(elevation*pi)/180; 
[xf,yf,zf|=sph2cart(vector_azimuth,vector_elevation,r); 

range=sqrt((position_x-xf).^2.+(position_y-yf).^2.+(position_z-zf).^2); 
% approximation 
angle=atan(range ./r); 

% matrixpos 

matrixpos=[r;phi;theta;angle]; 

% save to file, matrixpos 

s=sprintf('size%d. mat',i); 
save(s, 'matrixpos'); 


end 


%■ 


% add to E-field, drop and rain rate comparison 
sum_E=[sum_E,total_E]; 
weighted_E=[weighted_E,new_E]; 

end 


%■ 
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% - 

% Calculate RCS and doppler 


% calculate RCS, Z and Power return for every pulse including error 
running_Z=(Norm*abs(sum_E).^2)*wavelength^4/pi^5/abs(kw)^2/volume_resolution/... 
le-18+error_Z; 

power_retum=Norm*abs(weighted_E).^2; % [W] 

% calculate average Power and Z 
Z_avg=sum(mnning_Z)/p; 

Z_avg_plot=ones( 1 ,p) *Z_avg; 

power_retum_avg=sum(Norm*abs(weighted_E).^2)/p; % [W] 

power_retum_avg_plot=ones(l,p)*power_retum_avg; 

% calculate doppler frequency 
phase_vector_u=unwrap(phase_vector); 

doppler_fq=([phase_vector_u(2:2:2*p)-phase_vector_u(l:2:2*p-l)])/2/pi/PRT_l; 
number_vector=l :p; 


% - 

% Calculate references 
% Integral form of Z 
ii=0:100/1000:100; 

zz=ii.^6.*8000.*exp(-4.1.*rain_rate *(-.21).*ii); 
real_Z=ones( 1 ,number_pulses)* 10*log 10(trapz(ii,zz)); 

% Adding the error to the weighted E estimations 

error_power=error_rcs*P_t*G''2*wavelength^2/pi^2*(2*angle_3db*pi/l 80)^2*... 
(c*pulsewidth*le-6)/1024/log(2)/pi^2/(distance*le3)^2/volume_resolution; 

% Evaluate average Power return to make Z estimation including error 
Z_estimation=power_retum_avg/P_t/G^2/wavelength^2*... 

pi^2/(2*angle_3db*pi/180)^2/(c*pulsewidth*le-6)*1024*log(2)*(distance*le3)^2/... 
((kw''2)*(pi^5)/(wavelength^4))/le-18; 

Z_estimation_plot=ones( 1 ,number_pulses)* 10*log 10(Z_estimation+error_Z); 


% - 

% Part V, Plot of results 
%- 


figure(l); 

subplot(2,l,l); 

plot(number_vector(l :p),doppler_fq(l :p),'-x',number_vector(l :p),ones(l,p)*... 
mean(doppler_fq),r'); 

title(['Doppler frequency for all pulse pairs for PRF= ',num2sfr(l/PRT_l),' Hz']); 

xlabel('Pulse pair'); 

ylabel('f_d'); 

legend('Doppler frequency per pulse pair','Average Doppler estimation') 
grid on 

subplot(2,l,2) 
hist(doppler_fq, 15) 
hold on 

[ff,xx] = ksdensity(doppler_fq); 
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plot(xx,ff/max(ff)*max(hist(doppler_fq, 15))) 
title('Doppler Frequency Flistogram with Gaussian Fit'); 
ylabel('Frequency'); 
xlabel('Doppler Frequency [Fiz]'); 


figure (2) 
subplot(2,l,l); 

% also converting to dBm 

plot(number_vector*PRT_2*1000,10*logl0(power_retum+error_power)+30,'-xb',... 

number_vector*PRT_2*1000,10*logl0(power_retum_avg_plot+error_power)+30,'—r'); 
title(['Power Return for all ',num2str(p),' pulses. PRF= ',num2str(l/PRT_2),' Flz']); 
xlabel('Sample time [ms]'); 
ylabel('Power [dBm]'); 

legend('Power return', 'Average Power return',0) 
grid on 


subplot(2,l,2); 

plot(number_vector, 10*log 10(abs(running_Z)),'-xb',number_vector, 10*... 
loglO(abs(Z_avg_plot)),'—r',number_vector,real_Z,'—k',... 
number_vector,Z_estimation_plot,'—b'); 
title('Estimated Z per pulse'); 
xlabel('Sample time'); 
ylabel('Z [dBZ]'); 

legend('Simulation Z, (instantaneous)', 'Estimated Z, (average return)',... 

'Z=\int D^6N(D)dD','Estimated Z, Zeroth moment',0) 
grid on 

figure(3) 

polarmatrix(l :2:119)=zeros(l,60); 
polarmatrix(2:2:120)=phase_vector_u(l :60); 
lengthmatrix(l :2:119)=zeros(l,60); 
lengthmabix(2:2:120)=ones( 1,60); 
subplot(5,6,l) 

polar(polarmatrix(l :2),lengthmatrix(l :2),'-*') 
hold on 

polar(polarmatrix(3:4),lengthmatrix(3:4),'—*r') 

title(['l^{st}. v_r= ',num2str(doppler_fq(l)*wavelength/2),' m/s']) 

hbne = findobj(gca,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,2) 

polar(polarmatrix(5:6),lengthmatrix(5:6),'-*') 
hold on 

polar(polarmatrix(7:8),lengthmatrix(7:8),'—*r') 

title(['2^ {nd}. v_r= ',num2str(doppler_fq(2)*wavelength/2),' m/s']) 

hbne = findobj(gca,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,3) 

polar(polarmatrix(9:10),lengthmatrix(9:10),'-*') 
hold on 

polar(polarmatrix(l 1:12),lengthmatrix(l 1:12),':*r') 

title(['3^{rd}. v_r= ',num2str(doppler_fq(3)*wavelength/2),' m/s']) 

hbne = findobj(gca,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,4) 
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polar(polarmatrix(l 3:14),lengthmatrix(l 3:14),'-*') 
hold on 

polar(polarmatrix(l 5:16),lengthmatrix(l 5:16),':*r') 

title(['4^ {th}. v_r= ',num2str(doppler_fq(4)*wavelength/2),' m/s']) 

hline = findobj(gca,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,5) 

polar(polarmatrix( 17:18), lengthmatrix( 17:18),'-*') 
hold on 

polar(polarmatrix( 19:20), lengthmatrix( 19:20),':*r') 

title(['5^{th}. v_r= ',num2sh-(doppler_fq(5)*wavelength/2),' m/s']) 

hline = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,6) 

polar(polarmatrix(21:22),lengthmatrix(21:22),'-*') 
hold on 

polar(polarmatrix(23:24),lengthmatrix(23:24),':*r') 

title(['6^{th}. v_r= ',num2sh-(doppler_fq(6)*wavelength/2),' m/s']) 

hline = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,7) 

polar(polarmatrix(25:26),lengthmatrix(25:26),'-*') 
hold on 

polar(polarmatrix(27:28),lengthmatrix(27:28),':*r') 

title(['7^{th}. v_r= ',num2sh-(doppler_fq(7)*wavelength/2),' m/s']) 

hline = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,8) 

polar(polarmatrix(29:30),lengthmatrix(29:30),'-*') 
hold on 

polar(polarmatrix(31:32),lengthmatrix(31:32),':*r') 

title(['8^{th}. v_r= ',num2sb(doppler_fq(8)*wavelength/2),' m/s']) 

Mine = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,9) 

polar(polarmatrix(33:34),lengthmatrix(33:34),'-*') 
hold on 

polar(polarmatrix(3 5:3 6),lengthmatrix(3 5:3 6),': *r') 

title(['9^{th}. v_r= ',num2sb'(doppler_fq(9)*wavelength/2),' m/s']) 

Mine = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,10) 

polar(polarmatrix(3 7:3 8),lengthmatrix(3 7:38),'-*') 
hold on 

polar(polarmatrix(3 9:40),lengthmatrix(3 9:40),': *r') 

title(['10''{th}. v_r= ',num2sh-(doppler_fq(10)*wavelength/2),' m/s']) 

Mine = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,ll) 

polar(polarmatrix(41:42),lengthmatrix(41:42),'-*') 
hold on 

polar(polarmatrix(43:44),lengthmatrix(43:44),':*r') 

title(['ll^{th}. v_r= ',num2sb-(doppler_fq(ll)*wavelength/2),' m/s']) 

Mine = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,12) 
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polar(polarmatrix(45:46),lengthmatrix(45:46),'-*') 
hold on 

polar(polarmatrix(47:48),lengthmatrix(47:48),':*r') 

title(['12^{th}. v_r= ',num2str(doppler_fq(12)*wavelength/2),' m/s']) 

hline = findobj(gca,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,13) 

polar(polarmatrix(49:50), lengthmatrix(49:50),'- *') 
hold on 

polar(polarmatrix(51:52),lengthmatrix(51:52),':*r') 

title(['13^{th}. v_r= ',num2sh-(doppler_fq(13)*wavelength/2),' m/s']) 

hline = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,14) 

polar(polarmatrix(53:54),lengthmatrix(53:54),'-*') 
hold on 

polar(polarmatrix(5 5:56), lengthmatrix(5 5:5 6),': *r') 

title(['14^{th}. v_r= ',num2sb(doppler_fq(14)*wavelength/2),' m/s']) 

hline = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,15) 

polar(polarmatrix(57:58),lengthmatrix(57:58),'-*') 
hold on 

polar(polarmatrix(59:60),lengthmatrix(59:60),': *r') 

title(['15^{th}. v_r= ',num2sb(doppler_fq(15)*wavelength/2),' m/s']) 

hline = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,16) 

polar(polarmatrix(61:62),lengthmatrix(61:62),'-*') 
hold on 

polar(polarmatrix(63:64),lengthmatrix(63:64),':*r') 

title(['16^{th}. v_r= ',num2sh-(doppler_fq(16)*wavelength/2),' m/s']) 

Mine = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,17) 

polar(polarmatrix(65:66),lengthmatrix(65:66),'-*') 
hold on 

polar(polarmatrix(67:68),lengthmatrix(67:68),':*r') 

title(['17^{th}. v_r= ',num2sb-(doppler_fq(17)*wavelength/2),' m/s']) 

Mine = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,18) 

polar(polarmatrix(69:70),lengthmatrix(69:70),'-*') 
hold on 

polar(polarmatrix(71:72),lengthmatrix(71:72),':*r') 

title(['18^{th}. v_r= ',num2sh-(doppler_fq(18)*wavelength/2),' m/s']) 

Mine = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,19) 

polar(polarmatrix(73:74),lengthmatrix(73:74),'-*') 
hold on 

polar(polarmatrix(75:76), lengthmatrix(75:76),': *r') 

title(['19^{th}. v_r= ',num2sb-(doppler_fq(19)*wavelength/2),' m/s']) 

Mine = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,20) 
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polar(polarmatrix(77:78),lengthmatrix(77:78),'-*') 
hold on 

polar(polarmatrix(79:80), lengthmatrix(79:80),*r') 

title(['20^ {th}. v_r= ',num2str(doppler_fq(20)*wavelength/2),' m/s']) 

hline = findobj(gca,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,21) 

polar(polarmatrix(81:82),lengthmatrix(81:82),'-*') 
hold on 

polar(polarmatrix(83:84),lengthmatrix(83:84),':*r') 

title(['21^{th}. v_r= ',num2sh-(doppler_fq(21)*wavelength/2),' m/s']) 

hline = findobj(gca,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,22) 

polar(polarmatrix(85:86),lengthmatrix(85:86),'-*') 
hold on 

polar(polarmatrix(87:88),lengthmatrix(87:88),':*r') 

title(['22^{th}. v_r= ',num2stt(doppler_fq(22)*wavelength/2),' m/s']) 

hline = findobj(gca,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,23) 

polar(polarmatrix(89:90), lengthmatrix(89:90),'-*') 
hold on 

polar(polarmatrix(91:92),lengthmatrix(91:92),':*r') 

title(['23^{th}. v_r= ',num2str(doppler_fq(23)*wavelength/2),' m/s']) 

hline = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,24) 

polar(polarmatrix(93:94),lengthmatrix(93:94),'-*') 
hold on 

polar(polarmatrix(95:96),lengthmatrix(95:96),':*r') 

title(['24^ {th}. v_r= ',num2str(doppler_fq(24)*wavelength/2),' m/s']) 

Mine = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,25) 

polar(polarmatrix(97:98),lengthmatrix(97:98),'-*') 
hold on 

polar(polarmatrix(99:100),lengthmatrix(99:100),':*r') 

title(['25^{th}. v_r= ',num2sh-(doppler_fq(25)*wavelength/2),' m/s']) 

Mine = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,26) 

polar(polarmatrix(101:102),lengthmatrix(101:102),'-*') 
hold on 

polar(polarmatrix(103:104),lengthmatrix(103:104),':*r') 

title(['26^{th}. v_r= ',num2sh-(doppler_fq(26)*wavelength/2),' m/s']) 

Mine = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,27) 

polar(polarmatrix( 105:106),lengthmatrix( 105:106),'-*') 
hold on 

polar(polarmatrix( 107:108),lengthmatrix( 107:108),': *r') 

title(['27^ {th}. v_r= ',num2sh-(doppler_fq(27)*wavelength/2),' m/s']) 

Mine = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,28) 
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polar(polarmatrix( 109:110),lengthmatrix( 109:110),'-*') 
hold on 

polar(polan'natrix(l 11:112),lengthmatrix(l 11:112),':*r') 

title(['28''{th}. v_r= ',num2str(doppler_fq(28)*wavelength/2),' m/s']) 

hline = findobj(gca,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,29) 

polar(polan'natrix(l 13:114),lengthmatrix(l 13:114),'-*') 
hold on 

polar(polamiatrix(l 15:116),lengthmatrix(l 15:116),':*r') 

title(['29'''{th}. v_r= ',num2str(doppler_fq(29)*wavelength/2),' m/s']) 

hline = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,30) 

polar(polarmatrix(l 17:118),lengthmatrix(l 17:118),'-*') 
hold on 

polar(polarmatrix( 119:120),lengthmatrix( 119:120),': *r') 
title(['30''{th}. v_r= ',num2str(doppler_fq(30)*wavelength/2),' m/s']) 
hline = findobj(gea,'Type','line'); 
set(hline,'LineWidth',3) 

mean(doppler_fq) 

std(doppler_fq) 

3) WR_fc|_step 

% Weather Radar Signal Simulator, developed version 
% Frequeney step 
% Ulf Sehroder 
% August, 2005 

% - 


% Main program for The NPS/Marielle Gosset Weather Radar Signal Simulator. 


% - 

% DESCRIPTION 
% Part I 

% From Radar Data: 

% - Generate Beam Resolution eell approximation 

% - Generate Cube adding a margine to the Beam Resolution eell to allow for 
% the rain drops to fall trough the through beam over time. 

% 

% Part II 

% From Rain Parameters: 

% Generating the fundamental T Matrix eoneisting of 
% - Diameter Veetor and Number of drops in eaeh diameter interval 
% using the Marshall-Palmer exponential approximation. 

% - the total number of drops for eaeh size per m*3 
% The baekseatter RCSs for eaeh Diameter size is put in a separate RCS 
% Matrix to aeeount for the differenees due to frequeney 

% Part III 

% Randomly plaee all drops in eube 
% Part IV 
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% For every pulse the eoherent Eleetrie field return is ealeulated. Between 
% pulses all drops are moved. 

% Part V 

% Caleulate and plot results 

% - 

elear all 
elose all 
ele 

% Input Data 

e=3e8; % Speed of light [m/s] 

% Radar Data 

P_t= 10e3; % Transmitted power 

G_dB=20; % Maximum gain in dB 

G=10^(G_dB/10); % Maximum gain 

angle_3db=0.20; % 3 dB Beam angle [degrees], (half half power beamwidth) 

pulsewidth=0.20; % Pulsewidth [mieroseeond] 

elevation=3; % Elevation angle [degrees] 

margine=5; % margin for resolution eell[m] 

azimuth=0; % Azimuth angle [rad] 

number_pulses=40; % Number of pulses 

PRT_l=l/5000; % Doppler Pulse Repetition Time [seeonds] 

PRT_2=l/200; % Avg Power Pulse Repetition Time [seeonds] 

distanee=0.5; % Range from the radar to the resolution eell [Km] 
wavelength=0.1; % Wavelength [m] 

frequeney=c/wavelength; % Frequeney [Hz] 

% - 

%Introdueed in stepped version 
fq_step=le6; % Frequeney step in [Hz] 

delay=50e-6; % Delay for ehanging fq in [s] 

n_steps=3; % Number of Frequeney steps 

% - 

% Rain parameters 

diam_low=0.05; % Lowest diameter of Preeipitation [mm] 

diam_hm=l; % Divider between low and mid [mm] 

diam_mid=5; % Divider between mid and high [mm] 

diam_high=15; % Highest diameter of Preeipitation [mm] 
resolutionl=0.08; % Lower diameter resolution [mm] 
resolution2=0.3; % Higher diameter resolution [mm] 
resolution3=0.8; % Highest diameter resolution [mm] 
sort_approx=l; % RCS approximation model; [l]=Rayleigh 
kw=0.93; % abs(K)*2 approx for water 

rain_rate=50; % Rain Rate [mm/h] 

rain_anglel=0; % degrees Rain Angle (Phi) 

rain_angle2=0; % degrees Rain Angle (Theta) 
wind_veetor=[-10;0;0]; % Wind speed veetor (x,y,z) [m/s] 
spread=0; % Wind speed spread [m/s], (all direetions) 

spread_2=4; % Drop veloeity spread [m/s], 

Norm=100; % Normalization eoeffieient 

% Direetional eosines 
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u=sin(pi/2-elevatioii*pi/180)*cos(azimuth*pi/180); 
v=sin(pi/2-elevatioii*pi/l 80)*sin(azimuth*pi/l 80); 
w=cos(pi/2-elevatioii*pi/180); 


% - 

% Part I, Generate the resolution volume 

% - 

% Approximation of the resolution eell of the pulse eireular box, 

% for any azimuth angle 

% - 

% INPUT 
% - distanee 
% - angle_3dB 
% - pulsewidth 
% - elevation 
% - margine 
% - azimuth 

% - 

% OUTPUT 

% - volume resolution: eomputed volume of the resolution eell 
% - volume box: volume of the box where the drops are put 
% - box: the 8 reetangular eoordinates of the box 
% - eoord vol res: the eoordinates of the resolution eell where the drops 
% are taken. 

% - 


% - 

% Creating the eireular beam resolution eell 
range_res=(e*pulsewidth* le-6)/2; 

% Range resolution (et/2) [m] 

% small angle approximation allows for following 
angle_rad=(angle_3 db *pi)/180; 

% 3 dB beamwidth [rad] 
radius l=(tan(angle_rad)*(distanee* 1000)); 

% Radius of beam at R=distanee 
radius2=(tan(angle_rad)*((distanee* 1000)+range_res)); ... 

% Radius of beam at R=distanee + range resolution 

surfaee 1 =pi *radius U2; 
surfaee2=pi *radius2*2; 

% volume of eone = base surfaee times height for all three 
volumeeone 1 =(surfaee 1 *(distanee* 1000))/3; 
volume_eone2=(surfaee2*((distanee*1000)+range_res))/3; 
volume_resolution=(volume_eone2-volume_eonel); % Beam Cone Volume 

% eomputation of the resolution eell eoordinates, spherieal eoordinates 
% elevation in radians 

elevation_rad=(elevation*pi)/180; % Elevation [rad] 

azimuth_rad=(azimuth*pi)/180; % Azimuth [rad] 

% Matrix used when eomparing wheather a drop is inside or outside 
% beamwidth 

eoord_vol_res=[distanee* 1000,(distanee* 1000)+range_res;elevation_rad-angle_rad,... 
elevation_rad+angle_rad;azimuth_rad-angle_rad,azimuth_rad+angle_rad]; 
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%■ 


% - 

% Creating the box 

% eomputation of the box and its eoordinates 

eoinl=[azimuth_rad+angle_rad,elevation_rad+angle_rad,distanee* 1000]; 
eoin2=[azimuth_rad-angle_rad,elevation_rad+angle_rad,distanee* 1000]; 
eoin3=[azimuth_rad+angle_rad,elevation_rad-angle_rad,distanee *1000]; 
eoin4=[azimuth_rad-angle_rad,elevation_rad-angle_rad,distanee*1000]; 
eoin5=[azimuth_rad+angle_rad,elevation_rad+angle_rad,distanee* 1000+range_res]; 
eoin6=[azimuth_rad-angle_rad,elevation_rad+angle_rad,distanee*1000+range_res]; 
eoin7=[azimuth_rad+angle_rad,elevation_rad-angle_rad,distanee*1000+range_res]; 
eoin8=[azimuth_rad-angle_rad,elevation_rad-angle_rad,distanee* 1000+range_res]; 

[xl,yl,zl]=sph2eart(eoinl(l),eoinl(2),eoinl(3)); 

[x2,y2,z2]=sph2eart(eoin2(l),eoin2(2),eoin2(3)); 

[x3 ,y3 ,z3]=sph2eart(eoin3 (l),eoin3 (2),eoin3 (3)); 
[x4,y4,z4]=sph2eart(eoin4(l),eoin4(2),eoin4(3)); 
[x5,y5,z5]=sph2eart(eoin5(l),eoin5(2),eoin5(3)); 
[x6,y6,z6]=sph2eart(eoin6(l),eoin6(2),eoin6(3)); 

[x7,y7,z7]=sph2eart(eoin7 (l),eoin7 (2),eoin7 (3)); 
[x8,y8,z8]=sph2eart(eoin8(l),eoin8(2),eoin8(3)); 

xmax=max([xl,x2,x3,x4,x5,x6,x7,x8]); 

xmin=min([xl,x2,x3,x4,x5,x6,x7,x8]); 

ymax=max([yl,y2,y3,y4,y5,y6,y7,y8]); 

ymin=min([y 1 ,y2,y3 ,y4,y5,y6,y7,y8]); 

zmax=max([zl,z2,z3,z4,z5,z6,z7,z8]); 

zmin=min([zl,z2,z3,z4,z5,z6,z7,z8]); 

% the box with its margin 

eoin_l=[xmin-margine,ymin-margine,zmin-margine]; 

eoin_2=[xmin-margine,ymax+margine,zmin-margine]; 

eoin_3=[xmax+margine,ymax+margine,zmin-margine]; 

eoin_4=[xmax+margine,ymin-margine,zmin-margine]; 

eoin_5=[xmin-margine,ymin-margine,zmax+margine]; 

eoin_6=[xmin-margine,ymax+margine,zmax+margine]; 

eoin_7=[xmax+margine,ymax+margine,zmax+margine]; 

eoin_8=[xmax+margine,ymin-margine,zmax+margine]; 

box=[eoin r,eoin_2',eoin_3',eoin_4',eoin_5',eoin_6',eoin_7',eoin_8']; 

% ealeulate the volume of the box 

volume_box=abs((xmax+margine-(xmin-margine))*(ymax+margine-(ymin-margine))... 
* (zmax+margine-(zmin-margine))); 


eube=[eoinr,eoin_2',eoin_3',eoin_4',eoinr,eoin_5',eoin_6',eoin_2',... 

eoin_6',eoin_7',eoin_3',eoin_7',eoin_8',eoin_4',eoin_8',eoin_5']; 

% - 

% - 

% Part II, Generate the fundamental T Matrix 

% - 

% This part produees a matrix of all drops that eontribute to the total 
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% reflectivity of the resolution cell 

% - 

% First all sizes are created, with two vectors of different size resolution. 

% The size vectors represents all diameters. Next the backscatter RCS for each 
% size is computed. 

% Using the Marshall-Palmer exponential approximation, the number of drops of 
% each diameter is computed. 

% - 

% INPUT 
% - diam low 
% - diam lim 
% - diam high 
% - resolution! 

% - resolution2 
% - sort approx 
% - wavelength 
% - Kw 
% - rain_rate 

% - 

% OUTPUT 
% - number drops 
% Matrix T: 

% - first row all sizes 
% - second row all the backscatter RCSs 
% - third row the total number of drops for each size per m^3 
% the columns starting from the smallest size 

% - 

% Create dropsize vector 

% Precipitaion Diameter vector (lowest<D<divider) 
vector_size_ 1 =diam_low resolution 1 :diam_lim; 

% Precipitaion Diameter vector (divider<D<mid) 
vector_size_2=diam_lim:resolution2:diam_mid; 

% Precipitaion Diameter vector (mid<D<highest) 
vector_size_3=diam_mid:resolution3:diam_high; 

% Precipitaion Diameter Vector 

vector_size=[vector_size_l,vector_size_2,vector_size_3]; 


% - 

% Create RCS matrix for all Rain drops of Precipitation Vector and for all 
% frequencies. 

fq_rr=ff equency: fq_step: frequency+(n_steps-1 )* fq_step; 
rcs(l:n_steps,:)=((kw^2)*(pU5)/(cM)*fq_rr'.^4)*((vector_size*le-3).^6); 
%- 


% - 

% The Marshall Palmer approximation is used to create Precipitaion Diameter 
% Distribution Vector (Drop Size Distribution N(D) Vector) 

% The number of drops are scaled down by a factor (Norm) to improve the 
% speed of the simulator 

% - 

lambda=4.1 * (rain_rate^(-0.21)); % [mm-1 ] 

number=8000*exp(-lambda*vector_size)/Norm; 
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vectorres 1 =ones( 1 ,length(vector_size_l ))*resolutionl; 
vector_res2=ones(l,length(vector_size_2))*resolution2; 
vector_res3=ones(l,length(vector_size_3))*resolution3; 

% delta(Diameter) dD 

vector_res=[vector_resl,vector_res2,vector_res3]; 

% the approximate "true" number of drops N(D)dD vector 
number=number.*vector_res; 

% - 


%■ 


% create matrix T without res 
T=[vector_size;number]; 


%■ 


% - 

% Change matrix 'T' to matrix 'matrix T Round' 

% - 

% This part computes the total number of drops by utilizing the resolution 
% cell and the matrix T Round 

% The only difference between the matrix T Round and T, is the fact that 
% matrix T Round has the total number of drops on the third row and they 
% are rounded of to nearest integer value. 

% - 


matrix_T_Round=T; 

matrix_T_Round(2, :)=round(matrix_T_Round(2, :)*volume_box); 

% Keep track of errors due to round off 
diff=T (2,: )* volume_box-matrix_T_Round(2,:); 
error_rcs=diff*rcs'*Norm/volume_box*volume_resolution; 
error_Z=error_rcs/volume_resolution*c^4./fci_rr.^4/pi^5/abs(kw)^2/le-18; 


% - 

% Part III, Initial positioning of the drops within the volume 

% - 

% This part generates the initial position of the drops in the volume of 
% the box. 

% The positions are recorded in the file "sizexx.maf taking into account 
% the azimuth angle 

% - 

% INPUT 
% - box 

% - matrix T Round 
% - elevation 
% - azimuth 

% - 

% OUTPUT 

% All fdes including the sizes. Each contains the positions of all drops 
% for a given diameter. The data recoreded as follows: 

% - radial distance r 
% - elevation phi 
% - azimuth theta 

% - angle with respect to beam center. This parameter will be 

119 













% used to identify drops whieh are in the beam. 
%- 


% Randomly put the ealeulated number of drops of eaeh diameter in the box 
number_size=length(matrix_T_Round( 1,:; 

for i=l inumber size 

% seleet number of drops of eertain Diameter 
drops=matrix_T_Round(2,i); 

% ereate as many positions on the x axis as # drops per size 
elements_x=rand( 1,drops); 

% put them in the range of interest 
% Random plaeement x for every drop 

position_x=box( l,l)+(box(l,3)-box(l,l))* elementsx; 

% ereate as many positionson the y axis as # drops per size 
elements_y=rand( 1,drops); 

% put them in the range of interest 
% Random plaeement y for every drop 

position_y=box(2,l)+(box(2,2)-box(2,l))*elements_y; 

% ereate as many positions on the z axis as # drops per size 
elements_z=rand( 1,drops); 

% put them in the range of interest 
% Random plaeement z for every drop 

position_z=box(3,1 )+(box(3,5)-box(3,1 ))*elements_z; 

% matrix of results 

% Creating a drop Matrix using spherieal eoordinates 

[theta,phi,r]=eart2sph(position_x,position_y,position_z); 

% Caleulating the range, phi, theta and angle to every drop with 
% referenee to phi O and theta O. Used to estimate weighted power return. 
veetor_azimuth=ones( 1 ,drops) * ((azimuth *pi)/180); 
veetor_elevation=ones(l,drops)*((elevation*pi)/180); 
[xf,yf,zf|=sph2eart(veetor_azimuth,veetor_elevation,r); 
range=sqrt((position_x-xf). ^2 .+(position_y-y f) .^2. +(position_z-zf) .^2); 

% approximation 
angle=atan(range./r); 

% Saving the drop positions in size%d.mat 
% matrixpos 

matrixpos=[r;phi;theta;angle]; 
save(sprintf('size%d.mat',i),'matrixpos'); 
end 

% - 
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% - 

% Part IV, Calculate Coherent Electric field return. 
%- 


% For every pulse the coherent Electric field return is calculatet. Between 
% pulses all drops are moved. 


% 


% Run the simulation for specified number of pulses 
%- 


% Zero all vectors and matrices that are going to be used 

phase_vector=zeros(l,2*number_pulses); 

imaginary_retum=zeros(l,2*number_pulses); 

sum_E=zeros(n_steps,ceil(number_pulses/n_steps)); 

weighted_E=zeros(n_steps,ceil(number_pulses/n_steps)); 

running_Z=[]; 

% Initiate reference variables 

p=0; % Used to separate pulses 

hh=0; % Used to separate avg power samples 

while p<number_pulses 

step countei^O; % Used to separate frequencies 

hh=hh+l; 

for fc|_new=frequency: fq_step: frequency+(n_steps-1) * fq_step 

k=2*pi*fq_new/c; 

p=p+l; 

step_counter=step_counter+1; 

disp(sprintf('calculate the power-return of pulse %d\n',p)); 


% - 

% computes the (weighted) I and Q return 

% - 

% computes the total number of drops 

% - 

% INPUT 
% - Matrix T 
% - wavelength 
% - angle SdB 
% - distance 
% - coord vol res 
% - number_pulses 

% - 

% OUTPUT 

% - total E: the received complex Electric field) 

% - power: the power 
% - voltage: the voltage 

% - drops beam: the number of drops in the beam 

% - 

% Only the drops which are in the resolution cell are selected for the 
% computation 

% - 
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number_size=length(matrix_T_Round(l,:)); 

total_E=0; 

total_phase=0; 

new_E=0; 

for i=l :number_size 

% opening the file of registered data matrixpos. 
load(sprintf('size%d.mat',i)); 

% removing drops that are outside the resolution cell 
number_drops_per_size=length(matrixpos( 1,:; 
dropmatrix=[]; 

% Filling the drop matrix with all drops inside the beam 

% resolution volume 

for d=l :number_drops_per_size 

if ((matrixpos(4,d)<=(angle_3db*pi/180))&(matrixpos(l,d)>=... 
coord_vol_res( 1,1 ))&(matrixpos(l ,d)<=coord_vol_res( 1,2))); 
dropmatrix=[dropmatrix,matrixpos(: ,d)]; 
end 
end 


% - 

% Convert to cartesian coordinates 
if isempty(dropmatrix)==0 

[position_x,position_y,position_z]=sph2cart(dropmatrix(3, 

dropmatrix(2,:),dropmatrix(l,:)); 

% Create comparison vectors 

vector_azimuth=ones( 1, length(dropmatrix( 1,:))) * ((azimuth*pi)/180); 
vector_elevation=ones( 1 ,length(dropmatrix( 1,: )))* ((elevation*pi)/180); 


% Create a phase vector for the drops of current Diameter 
r_phase=position_x.*u+position_y.*v+position_z.*w; 
phase=exp((2 *j *k). *r_phase); 

% Sum the contributions to the E-field 

total_phase=total_phase+sum(phase); 

total_E=total_E+sum(sqrt(rcs(step_counter,i)).*phase); 


% - 

% ceate weighted power return, concidering range (r), theta, phi 
% with reference to boresight and radar parameters 
E_weighted=sqrt(P_t)*G*c/fq_rr(step_counter)/(sqrt(4*pi)^3)./... 
dropmatrix(l,:) *2.*(exp(-4*log(2).*((dropmatrix(3,:)-... 
vector_azimuth).^2/(2*angle_3db*pi/l 80)^2+((dropmatrix(2,:)... 
-vector_elevation) .^2/(2*angle_3db *pi/180)^2)))); 
new_E=new_E+sum(E_weighted. * (sqrt(rcs(step_counter,i)). *phase)); 

end 


%■ 
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phase_vector(2*p-1 )=ANGLE(total_phase); 
imaginary _retum(2 *p-1 )=total_phase; 

% add to E-field and drop- and rate-vector 

sum_E(step_counter,hh)=total_E; 

weighted_E(step_counter,hh)=new_E; 


%■ 


% - 

% This part simulates the movement of the drops and gives new positions. 
% First using the Doppler PRT (PRT l) 

% - 

% INPUT 

% - matrix T Round 
% - PRT 
% - azimuth 
% - elevation 
% - rain_anglel 
% - rain_angle2 
% - diam lim 
% - resolution! 

% - resolution2 

% - wind vector 

% - 

% OUTPUTS 

% Output files containing sizes with new positions 
% structured the same way as in previous 

% there is also an output fde for speed (one for each size) for the pulse 
% pair. 

% - 


number_size=length(matrix_T_Round(l,:)); 

% adding a Gaussian random term to make the wind speed spread 
if wind_vector(l)~=0 

new_wind_vector( 1 )=wind_vector( 1) *normmd( 1 ,abs(spread/max(wind_vector( 1)))); 
else 

new_wind_vector( 1 )=0; 
end 

if wind_vector(2)~=0 

new_wind_vector(2)=wind_vector(2)*normmd(l,abs(spread/max(wind_vector(2)))); 

else 

new_wind_vector(2)=0; 
end 

if wind_vector(3)~=0 

new_wind_vector(3 )=wind_vector(3) *normmd( 1 ,abs(spread/max(wind_vector(3)))); 
else 

new_wind_vector(3 )=0; 
end 

for i=l:number_size 

% opening of data file registered in matrixpos. 
load(sprintf('size%d.mat',i)); 
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% total number of drops 
drops=matrix_T_Round(2, i); 

% transforming to rectangular coordinates 
[x,y,z]=sph2cart(matrixpos(3,:),matrixpos(2,:),matrixpos(l,:)); 

diameter=matrix_T_Round( 1 ,i); 


% - 

% Drop fall speed depends on drop diameter in the z axis, 

% Atlas-Ulbrich approximation is used 

% Approximation is valid in the diameter range 5*e-4, 5*e-3 
%- 


wtmax=3 86.6* ((matrix_T_Round( l,i) *0.001)^0.67); 

velocity z=(wtmax*(ones(l,drops)))*cos((rain_anglel*pi)/l 80); 
velocityx=(wtmax*(ones(l,drops)))*sin((rain_anglel*pi)/180)*... 
cos((rain_angle2*pi)/l 80); 

velocityy=(wtmax*(ones(l,drops)))*sin((rain_anglel*pi)/180)*... 
sin((rain_angle2*pi)/l 80); 


% - 

% - 

% % adding another Gaussian random term to make the wind speed 

% % spread between different drops sizes 

% if new_wind_vector(l)~=0 

% newest_wind_vector( 1 )=new_wind_vector( 1) *normmd( 1 ,abs(spread_2/... 

% max(wind_vector(l)))); 

% else 

% newest_wind_vector(l)=0; 

% end 

% if new_wind_vector(2)~=0 

% newest_wind_vector(2)=new_wind_vector(2) *normmd( 1 ,abs(spread_2/... 

% max(wind_vector(2)))); 

% else 

% newest_wind_vector(2)=0; 

% end 

% if new_wind_vector(3)~=0 

% newest_wind_vector(3 )=new_wind_vector(3) *normmd( 1 ,abs(spread_2/... 

% max(wind_vector(3)))); 

% else 

% newest_wind_vector(3)=0; 

% end 

% 

% % If spread within pulse 

% velocityxt=newest_wind_vector(l)-velocityx; 

% velocityyt=newest_wind_vector(2)-velocityy; 

% velocityzt=newest_wind_vector(3)-velocityz; 

% 

% % speed registration 

% m_velocity=[velocityxt;velocityyt;velocityzt]; 

% save(sprintf('velocity%d.mat',i),'m_velocity'); 


%■ 
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% - 

% If spread for all drops regardless of size 
newest_wind_veetor=[]; 


% adding Gaussian random term to make the wind speed 
% spread between drops 
if new_wind_veetor(l)~=0 

newest_wind_veetor( 1, :)=new_wind_veetor( 1) *... 

normmd(l,abs(spread_2/max(new_wind_veetor(l))),l,drops); 

else 

newest_wind_veetor( 1, :)=0; 
end 

if new_wind_veetor(2)~=0 

newest_wind_veetor(2, :)=new_wind_veetor(2)*... 

normmd(l,abs(spread_2/max(new_new_wind_veetor(2))),l,drops); 

else 

newest_wind_veetor(2, :)=0; 
end 

if new_wind_veetor(3)~=0 
newest_wind_veetor(3, :)=new_wind_veetor(3) *... 

normmd(l ,abs(spread_2/max(new_wind_veetor(3))), 1,drops); 

else 

newest_wind_veetor(3, :)=0; 
end 

% If spread within pulse 
veloeityxt=newest_wind_veetor( 1, :)-velocityx; 
veloeityyt=newest_wind_veetor(2,:)-velocityy; 
veloeityzt=newest_wind_veetor(3,:)-veloeityz; 


%■ 

%■ 


deplaeementx=veloeityxt*PRT_l; 
deplaeementy=veloeityyt*PRT_l; 
deplaeementz=veloeityzt*PRT_l; 

position_x=x+deplaeementx; 
position_y=y+deplaeementy; 
po sition_z=z+dep lae ementz; 

% transformation into spherieal eoordinates 

[theta,phi,r]=eart2sph(position_x,position_y,position_z); 

% computation of new angles 
vector_azimuth=ones(l,drops)*(azimuth*pi)/180; 
vector_elevation=ones( 1, drops) * (elevation*pi)/180; 
[xf,yf,zf|=sph2cart(vector_azimuth,vector_elevation,r); 

range=sqrt((position_x-xf). ''2 .+(position_y-y 1) .^2 .+(position_z-zf). ^2); 
% approximation 
angle=atan(range. /r); 

% matrixpos 

matrixpos=[r;phi;theta;angle]; 
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% save to file, matrixpos 

s=spriiitf('size%d. mat',i); 
save(s, 'matrixpos'); 


end 


% - 

% Caleulate phase for pulse pair 
total_phase=0; 

for i=l mumber size 

% opening the file of registered data matrixpos. 
load(sprintf('size%d.mat',i)); 

% removing drops that are outside the resolution eell 
number_drops_per_size=length(matrixpos( 1,:; 
dropmatrix=[]; 

% Filling the drop matrix with all drops inside the beam 

% resolution volume 

for d=l :number_drops_per_size 

if ((matrixpos(4,d)<=(angle_3db*pi/180))&(matrixpos(l,d)>=... 
eoord_vol_res( 1,1 ))&(matrixpos(l ,d)<=coord_vol_res( 1,2))); 
dropmatrix=[dropmatrix,matrixpos(: ,d)]; 
end 
end 


% - 

% Convert to eartesian eoordinates 
if isempty(dropmatrix)==0 

[position_x,position_y,position_z]=sph2eart(dropmatrix(3, 

dropmatrix(2,:),dropmatrix(l,:)); 

% Create eomparison veetors 

veetor_azimuth=ones( 1 ,length(dropmatrix( 1,:))) * ((azimuth*pi)/180); 
veetor_elevation=ones(l,length(dropmatrix(l,:)))*((elevation*pi)/180); 


% Create a phase veetor for the drops of eurrent Diameter 
r_phase=position_x. *u+position_y. *v+position_z. *w; 
phase=exp((2 *j *k). *r_phase); 

% Sum the eontributions to the E-field 
total_phase=total_phase+sum(phase); 
end 
end 

phase_veetor(2*p)=ANGLE(total_phase); 
imaginary _retum(2*p)=total_phase; 


% - 

% This part simulates the movement of the drops and gives new positions. 
% Using the doppler PRT and the delay (PRT_2, delay) 
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%■ 


%■ 


number_size=length(matrix_T_Round(l,:)); 

% adding a Gaussian random term to make the wind speed spread 
if wind_vector(l)~=0 

new_wind_vector(l)=wind_vector(l)*normmd(l,abs(spread/... 
max(wind_vector( 1)))); 

else 

new_wind_vector( 1 )=0; 
end 

if wind_vector(2)~=0 

new_wind_vector(2)=wind_vector(2)*normmd(l,abs(spread/... 
max(wind_vector(2)))); 

else 

new_wind_vector(2)=0; 
end 

if wind_vector(3)~=0 

new_wind_vector(3)=wind_vector(3)*normmd(l,abs(spread/... 
max(wind_vector(3)))); 

else 

new_wind_vector(3 )=0; 
end 

for i=l:number_size 

% opening of data file registered in matrixpos. 
load(sprintf('size%d.maf,i)); 

% total number of drops 
drops=matrix_T_Round(2,i); 

% transforming to rectangular coordinates 
[x,y,z]=sph2cart(matrixpos(3,:),matrixpos(2,:),matrixpos(l,:)); 

diameter=matrix_T_Round( 1 ,i); 


% - 

% Drop fall speed depends on drop diameter in the z axis, 

% Atlas-Ulbrich approximation is used 

% Approximation is valid in the diameter range 5*e-4, 5*e-3 
%- 


wtmax=386.6*((matrix_T_Round(l,i)*0.001)^0.67); 

velocityz=(wtmax*(ones(l,drops)))*cos((rain_anglel*pi)/180); 

velocityx=(wtmax*(ones(l,drops)))*sin((rain_anglel*pi)/180)*... 

cos((rain_angle2*pi)/180); 

velocityy=(wtmax* (ones( 1, drops))) * sin((rain_anglel*pi)/180)*... 
sin((rain_angle2 *pi)/180); 


%■ 

%■ 
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% % adding another Gaussian random term to make the wind speed 

% % spread between different drops sizes 

% i f ne w_wind_veetor( 1 )~=0 

% newest_wind_veetor( 1 )=new_wind_veetor( 1) *normmd( 1 ,abs(spread_2/... 

% max(wind_veetor(l)))); 

% else 

% newest_wind_veetor(l)=0; 

% end 

% if new_wind_vector(2)~=0 

% newest_wind_veetor(2)=new_wind_veetor(2) *normmd( 1 ,abs(spread_2/... 

% max(wind_veetor(2)))); 

% else 

% newest_wind_veetor(2)=0; 

% end 

% if new_wind_veetor(3)~=0 

% newest_wind_veetor(3 )=new_wind_veetor(3) *normmd( 1 ,abs(spread_2/... 

% max(wind_veetor(3)))); 

% else 

% newest_wind_veetor(3)=0; 

% end 

% 

% % If spread within pulse 

% veloeityxt=newest_wind_veetor( 1 )-veloeityx; 

% veloeityyt=newest_wind_veetor(2)-veloeityy; 

% veloeityzt=newest_wind_veetor(3)-veloeityz; 

% 

% - 

% - 

% If spread for all drops regardless of size 
newest_wind_veetor=[]; 

% adding Gaussian random term to make the wind speed 
% spread between drops 
if new_wind_veetor( 1 )~=0 

newest_wind_veetor( 1, :)=new_wind_veetor( 1) *... 

normmd(l,abs(spread_2/max(new_wind_veetor(l))),l,drops); 

else 

newest_wind_veetor( 1, :)=0; 
end 

if new_wind_veetor(2)~=0 

newest_wind_veetor(2,:)=new_wind_veetor(2)*... 

normmd(l,abs(spread_2/max(new_new_wind_veetor(2))),l,drops); 

else 

newe st_wind_veetor(2,: )=0; 
end 

if new_wind_veetor(3)~=0 

newest_wind_veetor(3, :)=new_wind_veetor(3)*... 

normmd(l ,abs(spread_2/max(new_wind_veetor(3))), 1,drops); 

else 

newest_wind_veetor(3, :)=0; 
end 
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% If spread within pulse 
veloeityxt=newest_wind_veetor(l,:)-velocityx; 
veloeityyt=newest_wind_veetor(2,:)-veloeityy; 
veloeityzt=newest_wind_veetor(3,:)-veloeityz; 


%■ 

%■ 


% speed registration 

m_veloeity=[veloeityxt;veloeityyt;veloeityzt]; 

save(sprintf('veloeity%d.mat',i),'m_veloeity'); 

deplaeementx=veloeityxt * (PRT_ 1 +delay); 
deplaeementy=veloeityyt* (PRT_ 1 +delay); 
deplaeementz=veloeityzt* (PRT_ 1 +delay); 

position_x=x+deplaeementx; 

position_y=y+deplaeementy; 

position_z=z+deplaeementz; 

% transformation into spherieal eoordinates 

[theta,phi,r]=eart2sph(position_x,position_y,position_z); 

% eomputation of new angles 
veetor_azimuth=ones( 1 ,drops) * (azimuth*pi)/180; 
veetor_elevation=ones( 1 ,drops)* (elevation*pi)/180; 
[xf,yf,zf|=sph2eart(veetor_azimuth,veetor_elevation,r); 

range=sqrt((position_x-xf) .^2 .+(position_y-y 1). ''2 .+(position_z-zf) .^2); 

% approximation 

angle=atan(range./r); 

% matrixpos 

matrixpos=[r;phi;theta;angle]; 

% save to file, matrixpos 

s=sprintf('size%d. mat',i); 
save(s, 'matrixpos'); 


end 


%■ 


% - 

% This part simulates the movement of the drops and gives new positions. 
% Using the Avg Power PRT redueed by the sample times for the previous 
% samples 

% - 


% - 

number_size=length(matrix_T_Round(l,:)); 

% adding a Gaussian random term to make the wind speed spread 
if wind_veetor(l)~=0 
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new_wind_vector( 1 )=wind_vector( 1 )*normmd(l ,abs(spread/... 
max(wind_vector( 1)))); 

else 

new_wind_veetor( 1 )=0; 
end 

if wind_veetor(2)~=0 

new_wind_veetor(2)=wind_veetor(2)*normmd(l,abs(spread/... 
max(wind_veetor(2)))); 

else 

new_wind_veetor(2)=0; 
end 

if wind_veetor(3)~=0 

new_wind_veetor(3 )=wind_veetor(3 )*normmd( 1 ,abs(spread/... 
max(wind_veetor(3)))); 

else 

new_wind_veetor(3 )=0; 
end 


for i=l :number_size 

% opening of data file registered in matrixpos. 
load(sprintf('size%d.mat',i)); 

% total number of drops 
drops=matrix_T_Round(2, i); 

% transforming to reetangular eoordinates 
[x,y,z]=sph2eart(matrixpos(3,:),matrixpos(2,:),matrixpos(l,:)); 

diameter=matrix_T_Round( 1 ,i); 


% - 

% Drop fall speed depends on drop diameter in the z axis, 

% Atlas-Ulbrieh approximation is used 

% Approximation is valid in the diameter range 5*e-4, 5*e-3 
%- 


wtmax=3 86.6* ((matrix_T_Round( l,i)*0.001)^0.67); 

veloeityz=(wtmax*(ones(l,drops)))*eos((rain_anglel*pi)/180); 
veloeityx=(wtmax*(ones(l,drops)))*sin((rain_anglel*pi)/180)*... 
eos((rain_angle2*pi)/l 80); 

veloeityy=(wtmax* (ones( 1,drops))) * sin((rain_anglel*pi)/180)*... 
sin((rain_angle2*pi)/l 80); 


% - 

% - 

% % adding another Gaussian random term to make the wind speed 

% % spread between different drops sizes 

% if new_wind_veetor(l)~=0 

% newest_wind_veetor( 1 )=new_wind_veetor( 1) *normmd( 1 ,abs(spread_2/... 

% max(wind_veetor(l)))); 

% else 

% newest_wind_veetor(l)=0; 

% end 
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% if new_wind_vector(2)~=0 

% newest_wind_vector(2)=new_wind_vector(2) *normmd( 1 ,abs(spread_2/. 

% max(wind_vector(2)))); 

% else 

% newest_wind_veetor(2)=0; 

% end 

% if new_wind_veetor(3)~=0 

% newest_wind_veetor(3 )=new_wind_veetor(3) *normmd( 1 ,abs(spread_2/. 

% max(wind_veetor(3)))); 

% else 

% newest_wind_veetor(3)=0; 

% end 

% 

% % If spread within pulse 

% veloeityxt=newest_wind_veetor(l)-veloeityx; 

% veloeityyt=newest_wind_veetor(2)-veloeityy; 

% veloeityzt=newest_wind_veetor(3)-veloeityz; 


% - 

% - 

% If spread for all drops regardless of size 
newest_wind_veetor=[]; 


% adding Gaussian random term to make the wind speed 
% spread between drops 
if new_wind_veetor(l)~=0 

newest_wind_veetor( 1, :)=new_wind_veetor( 1) *... 

normmd(l,abs(spread_2/max(new_wind_veetor(l))),l,drops); 

else 

newest_wind_veetor( 1, :)=0; 
end 

if new_wind_veetor(2)~=0 

newest_wind_veetor(2,:)=new_wind_veetor(2)*... 

normmd(l ,abs(spread_2/max(new_new_wind_veetor(2))), 1,drops); 

else 

newest_wind_veetor(2, :)=0; 
end 

if new_wind_veetor(3)~=0 
newest_wind_veetor(3, :)=new_wind_veetor(3) *... 

normmd(l ,abs(spread_2/max(new_wind_veetor(3))), 1,drops); 

else 

newest_wind_veetor(3, :)=0; 
end 

% If spread within pulse 
veloeityxt=newest_wind_veetor( 1, :)-veloeityx; 
veloeityyt=newest_wind_veetor(2,:)-veloeityy; 
veloeityzt=newest_wind_veetor(3,:)-veloeityz; 


%■ 

%■ 


% speed registration 

m_veloeity=[veloeityxt;veloeityyt;veloeityzt]; 
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save(sprintf('velocity%d.mat',i),'m_velocity'); 

deplacementx=velocityxt* (PRT_2-n_steps* (PRT_ 1 +delay)); 
deplacementy=velocityyt*(PRT_2-ii_steps*(PRT_l+delay)); 
deplacementz=velocityzt*(PRT_2-ii_steps*(PRT_l+delay)); 

positioii_x=x+deplacemeiitx; 

positioii_y=y+deplacementy; 

positioii_z=z+deplacementz; 

% transformation into spherical coordinates 

[theta,phi,r]=cart2sph(position_x,position_y,position_z); 

% computation of new angles 
vector_azimuth=ones(l,drops)*(azimuth*pi)/180; 
vector_elevation=ones( 1, drops) * (elevation*pi)/180; 
[xf,yf,zf|=sph2cart(vector_azimuth,vector_elevation,r); 

range=sqrt((position_x-xf). *2 .+(position_y-y f) .*2 .+(position_z-zf). ^2); 
% approximation 
angle=atan(range. /r); 

% matrixpos 

matrixpos=[r;phi;theta;angle]; 

% save to file, matrixpos 

s=sprintf('size%d. mat',i); 
save(s, 'matrixpos'); 
end 

end 

end 


%■ 


% - 

% Calculate RCS and doppler 
%- 


% - 

% calculate RCS, Z and Power return for every pulse including error 

% - 

for 0=1 :n_steps 

mnning_Z(o,:)=(Norm*abs(sum_E(o,:)).^2)*cM/fc|_rr(o)''4/pC5/... 

abs(kw)''2/volume_resolution/le-18+error_Z(o); 

end 

running_Z=reshape(running_Z,l,[]); % Transform into a vector 

% - 

% Adding the error to the weighted E estimations 
%- 


error_power=error_rcs*P_t*G''2*c''2./fq_rr.^2/pC2*(2*angle_3db*pi/l 80)^2*... 
(c*pulsewidth* 1 e-6)... 

/1024/log(2)/pi^2/(distance* 1 e3)^2; 
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for 0=1 :n_steps 

power_retum_matrix(o,:)=Norm*abs(weighted_E(o,:)).^2+error_power(o); 

end 

power_retum=reshape(power_retum_matrix,!,[]);% Transform into a vector 


% - 

% calculate average Power and Z 

% - 

Z_avg=sum(mnning_Z)/p; 

Z_avg_plot=ones( 1 ,p) *Z_avg; 

power_retum_avg=sum(sum(Norm*abs(weighted_E).^2,2)+error_power')/p; % [W] 
power_retum_avg_plot=ones(l,p)*power_retum_avg; 


% - 

% calculate doppler frequency 
% use unwrap to make sure no alias is presented 

% - 

phase_vector_u=imwrap(phase_vector); 

doppler_fq=([phase_vector_u(2:2:2*p)-phase_vector_u(l:2:2*p-l)])/2/pi/PRT_l; 
number_vector=l :p; 


% - 

% Calculate references 
% Integral form of Z 

% - 

ii=0:100/1000:100; 

zz=ii.''6.*8000.*exp(-4.1.*rain_rate.^(-.21).*ii); 
real_Z=ones( 1 ,p) * 10*log 10(trapz(ii,zz)); 


% - 

% Evaluate average Power return to make Z estimation including error 

% - 

% To make the estimate with the correct frequency, the power return Matrix 
% is used. For the plot the result must me summed and divided by the number 
% of samples used to average 

Z_estimation=power_retum_matrix'*(fq_rr.^2)'/P_t/G^2/c''2*... 
pi^2/(2*angle_3db*pi/180)^2/(c*pulsewidth*le-6)*1024*log(2)*(distance*le3)^2/... 
((kw^2)*(pi^5)/(wavelength^4))/le-18; 

Z_estimation_plot=ones(l,p)*10*logl0((sum(Z_estimation)+sum(error_Z))/p); 


% - 

% Part V, Plot of results 
%- 


figure(l); 

subplot(2,l,l); 

plot(number_vector(l :p),doppler_fq(l :p),'-x',number_vector(l :p),ones(l,p)*... 
mean(doppler_fq),r'); 

title(['Doppler frequency for all pulse pairs for PRF= ',num2sfr(l/PRT_l),' Hz']); 

xlabel('Pulse pair'); 

ylabel('f_d'); 
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legend('Doppler frequency per pulse pairVAverage Doppler estimation') 
grid on 


subplot(2,l,2) 

% To account for the differences in wavelengths a correct wavelength vector 
% is created and multiplied with the doppler frequency vector 
u=ones( 1 ,p/n_steps); 

lambida=[u*c/fq_rr(l);u*c/fq_rr(2);u*c/fq_rr(3)]; 
wave_l=reshape(lambida, 1,[]); 
v_r=doppler_fq. *wave_l/2; 

hist(v_r,15) 
hold on 

[ff,xx] = ksdensity(v_r); 

plot(xx,ff/max(ff)*max(hist(v_r,15))) 

title('Radial Velocity Histogram with Gaussian Fit'); 

xlabel('Radial velocity [m/s]'); 

ylabel('Frequency'); 

grid on 


figure(2) 

subplot(2,l,l); 

kl=0:p/n_steps-l; 

power_number_vector=[kl*PRT_2; kl*PRT_2+PRT_l+PRT_l+delay; kl*PRT_2+2*... 
(PRT_ 1 +PRT_ 1 +delay)]; 

power_number_vector=reshape(power_number_vector, 1,[]); 

text_string=['Power Return for all ',num2str(p),' pulses.',' Frequency stepped from 
num2str(fq_rr(l)/le9),' GHz to ',num2str(fq_rr(n_steps)/le9),' GHz with 
num2str(fq_step/le6),' MHz steps.']; 

% also converting to dBm 

plot(power_number_vector*1000,10*logl0(power_retum)+30,'-xb',... 

power_number_vector*1000,10*logl0(power_retum_avg_plot)+30,'—r'); 
title(text_string); 
xlabel('Sample time [ms]'); 
ylabel('Power [dBm]'); 

legend('Power return', 'Average Power retum',0) 
grid on 

subplot(2,l,2); 

plot(power_number_vector* 1000,10*log 10(abs(mnning_Z)),'-xb',power_number_vector*... 
1000,10*logl0(abs(Z_avg_plot)),'-.r',power_number_vector*1000,real_Z,'-k',... 
power_number_vector*1000,Z_estimation_plot,'—b'); 
title('Estimated Z per pulse'); 
xlabel('Sample time'); 
ylabel('Z [dBZ]'); 

legend('Simulation Z, (instantaneous)', 'Estimated Z, (average return)',... 

'Z=\int D^6N(D)dD','Estimated Z, Zeroth momenf,0) 
grid on 

figure (3) 

% This plot only works if the number of pulse pairs exceed 30 
polarmatrix(l :2:119)=zeros(l,60); 


134 



polarmatrix(2:2:120)=phase_vector_u(l :60); 
lengthmatrix(l :2:119)=zeros(l,60); 
lengthmatrix(2:2:120)=ones( 1,60); 
subplot(5,6,l) 

polar(polarmatrix(l :2),lengthmatrix(l :2),'-*') 
hold on 

polar(polarmatrix(3:4),lengthmatrix(3:4),'—*r') 

title(['l^{st}. v_r= ',num2str(doppler_fq(l)*c/fq_rr(l)/2),' m/s']) 

hline = findobj(gca,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,2) 

polar(polarmatrix(5:6),lengthmatrix(5:6),'-*') 
hold on 

polar(polarmatrix(7:8),lengthmatrix(7:8),'—*r') 

title(['2^ {nd}. v_r= ',num2str(doppler_fq(2)*e/fq_rr(2)/2),' m/s']) 

hline = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,3) 

polar(polarmatrix(9:10),lengthmatrix(9:10),'-*') 
hold on 

polar(polarmatrix(l 1:12),lengthmatrix(l 1:12),':*r') 

title(['3^{rd}. v_r= ',num2str(doppler_fq(3)*e/fq_rr(3)/2),' m/s']) 

hline = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,4) 

polar(polarmatrix( 13:14), lengthmatrix( 13:14),'-*') 
hold on 

polar(polarmatrix(l 5:16),lengthmatrix(l 5:16),':*r') 

title(['4^{th}. v_r= ',num2str(doppler_fq(4)*e/fq_rr(l)/2),' m/s']) 

hline = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,5) 

polar(polarmatrix( 17:18), lengthmatrix( 17:18),'-*') 
hold on 

polar(polarmatrix( 19:20), lengthmatrix( 19:20),':*r') 

title(['5^{th}. v_r= ',num2str(doppler_fq(5)*e/fq_rr(2)/2),' m/s']) 

hline = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,6) 

polar(polarmatrix(21:22),lengthmatrix(21:22),'-*') 
hold on 

polar(polarmatrix(23:24),lengthmatrix(23:24),':*r') 

title(['6^{th}. v_r= ',num2str(doppler_fq(6)*e/fq_rr(3)/2),' m/s']) 

hline = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,7) 

polar(polarmatrix(25:26),lengthmatrix(25:26),'-*') 
hold on 

polar(polarmatrix(27:28),lengthmatrix(27:28),':*r') 

title(['7^{th}. v_r= ',num2str(doppler_fq(7)*e/fq_rr(l)/2),' m/s']) 

hline = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,8) 

polar(polarmatrix(29:30),lengthmatrix(29:30),'-*') 
hold on 

polar(polarmatrix(31:32),lengthmatrix(31:32),':*r') 
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title(['8^{th}. v_r= ',num2str(doppler_fq(8)*c/fc|_rr(2)/2),' m/s']) 

hline = findobj(gca,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,9) 

polar(polarmatrix(33:34),lengthmatrix(33:34),'-*') 
hold on 

polar(polarmatrix(3 5:3 6), lengthmatrix(3 5:3 6),': *r') 

title(['9^{th}. v_r= ',num2str(doppler_fq(9)*e/fq_rr(3)/2),' m/s']) 

hline = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,10) 

polar(polarmatrix(3 7:3 8), lengthmatrix(3 7:38),'-*') 
hold on 

polar(polarmatrix(3 9:40),lengthmatrix(3 9:40),': *r') 

title(['10^{th}. v_r= ',num2str(doppler_fq(10)*e/fq_rr(l)/2),' m/s']) 

hline = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,ll) 

polar(polarmatrix(41:42),lengthmatrix(41:42),'-*') 
hold on 

polar(polarmatrix(43:44),lengthmatrix(43:44),':*r') 

title(['ll^{th}. v_r= ',num2str(doppler_fq(ll)*e/fq_rr(2)/2),' m/s']) 

hline = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,12) 

polar(polarmatrix(45:46),lengthmatrix(45:46),'-*') 
hold on 

polar(polarmatrix(47:48),lengthmatrix(47:48),':*r') 

title(['12^{th}. v_r= ',num2str(doppler_fq(12)*e/fq_rr(3)/2),' m/s']) 

hline = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,13) 

polar(polarmatrix(49:50), lengthmatrix(49:50),'- *') 
hold on 

polar(polarmatrix(51:52),lengthmatrix(51:52),':*r') 

title(['13^{th}. v_r= ',num2str(doppler_fq(13)*e/fq_rr(l)/2),' m/s']) 

hline = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,14) 

polar(polarmatrix(53:54),lengthmatrix(53:54),'-*') 
hold on 

polar(polarmatrix(5 5:56),lengthmatrix(5 5:5 6),': *r') 

title(['14^{th}. v_r= ',num2str(doppler_fq(14)*e/fq_rT(2)/2),' m/s']) 

hline = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,15) 

polar(polarmatrix(57:58),lengthmatrix(57:58),'-*') 
hold on 

polar(polarmatrix(59:60), lengthmatrix(59:60),': *r') 

title(['15^{th}. v_r= ',num2str(doppler_fq(15)*e/fq_rr(3)/2),' m/s']) 

hline = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,16) 

polar(polarmatrix(61:62),lengthmatrix(61:62),'-*') 
hold on 

polar(polarmatrix(63:64),lengthmatrix(63:64),':*r') 
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title(['16^{th}. v_r= ',num2str(doppler_fq(16)*c/fq_rT(l)/2),' m/s']) 

hline = findobj(gca,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,17) 

polar(polarmatrix(65:66),lengthmatrix(65:66),'-*') 
hold on 

polar(polarmatrix(67:68),lengthmatrix(67:68),':*r') 

title(['17^{th}. v_r= ',num2str(doppler_fq(17)*e/fc|_rr(2)/2),' m/s']) 

hline = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,18) 

polar(polarmatrix(69:70),lengthmatrix(69:70),'-*') 
hold on 

polar(polarmatrix(71:72),lengthmatrix(71:72),':*r') 

title(['18^{th}. v_r= ',num2str(doppler_fq(18)*e/fq_rr(3)/2),' m/s']) 

hline = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,19) 

polar(polarmatrix(73:74),lengthmatrix(73:74),'-*') 
hold on 

polar(polarmatrix(75:76), lengthmatrix(75:76),': *r') 

title(['19^{th}. v_r= ',num2str(doppler_fq(19)*e/fq_rr(l)/2),' m/s']) 

hline = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,20) 

polar(polarmatrix(77:78), lengthmatrix(77:78),'-*') 
hold on 

polar(polarmatrix(79:80), lengthmatrix(79:80),': *r') 

title(['20^ {th}. v_r= ',num2str(doppler_fq(20)*e/fq_rT(2)/2),' m/s']) 

hline = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,21) 

polar(polarmatrix(81:82),lengthmatrix(81:82),'-*') 
hold on 

polar(polarmatrix(83:84),lengthmatrix(83:84),':*r') 

title(['21^{th}. v_r= ',num2str(doppler_fq(21)*e/fq_rr(3)/2),' m/s']) 

hline = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,22) 

polar(polarmatrix(85:86),lengthmatrix(85:86),'-*') 
hold on 

polar(polarmatrix(87:88),lengthmatrix(87:88),':*r') 

title(['22^{th}. v_r= ',num2str(doppler_fq(22)*e/fq_rT(l)/2),' m/s']) 

hline = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,23) 

polar(polarmatrix(89:90),lengthmatrix(89:90),'- *') 
hold on 

polar(polarmatrix(91:92),lengthmatrix(91:92),':*r') 

title(['23^{th}. v_r= ',num2str(doppler_fq(23)*e/fq_rr(2)/2),' m/s']) 

hline = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,24) 

polar(polarmatrix(93:94),lengthmatrix(93:94),'-*') 
hold on 

polar(polarmatrix(95:96),lengthmatrix(95:96),':*r') 
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title(['24''{th}. v_r= ',num2str(doppler_fq(24)*c/fc|_rr(3)/2),' m/s']) 

hline = findobj(gca,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,25) 

polar(polarmatrix(97:98),lengthmatrix(97:98),'-*') 
hold on 

polar(polan'natrix(99:100),lengthmatrix(99:100),':*r') 

tide(['25''{th}. v_r= ',num2str(doppler_fq(25)*e/fq_rr(l)/2),' m/s']) 

hline = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,26) 

polar(polan'natrix(101:102),lengthmatrix(101:102),'-*') 
hold on 

polar(polarmatrix( 103:104),lengthmatrix( 103:104),': *r') 

tide(['26''{th}. v_r= ',num2str(doppler_fq(26)*e/fq_rr(2)/2),' m/s']) 

hline = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,27) 

polar(polarmatrix(105:106),lengthmatrix(105:106),'-*') 
hold on 

polar(polarmatrix( 107:108),lengthmatrix( 107:108),': *r') 

title(['27''{th}. v_r= ',num2str(doppler_fq(27)*e/fq_rr(3)/2),' m/s']) 

hline = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,28) 

polar(polannatrix( 109:110),lengthmatrix( 109:110),'-*') 
hold on 

polar(polan'natrix(l 11:112),lengthmatrix(l 11:112),':*r') 

title(['28''{th}. v_r= ',num2str(doppler_fq(28)*e/fq_rr(l)/2),' m/s']) 

hline = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,29) 

polar(polan'natrix(l 13:114),lengthmatrix(l 13:114),'-*') 
hold on 

polar(polamiatrix(l 15:116),lengthmatrix(l 15:116),':*r') 

title(['29''{th}. v_r= ',num2str(doppler_fq(29)*e/fq_rr(2)/2),' m/s']) 

hline = findobj(gea,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,30) 

polar(polarmatrix(l 17:118),lengthmatrix(l 17:118),'-*') 
hold on 

polar(polarmatrix( 119:120),lengthmatrix( 119:120),': *r') 
tide(['30''{th}. v_r= ',num2str(doppler_fq(30)*e/fq_rr(3)/2),' m/s']) 
hline = findobj(gea,'Type','line'); 
set(hline,'LineWidth',3) 


mean(doppler_fq) 

std(doppler_fq) 

4) WR_pulse_com 

% Weather Radar Signal Simulator, developed version 
% Pulse eompression 
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% Ulf Schroder 
% August, 2005 

% - 

% Main program for The NPS/Marielle Gosset Weather Radar Signal Simulator. 

% - 

% DESCRIPTION 
% Part I 

% From Radar Data: 

% - Generate Beam Resolution cell approximation 

% - Generate Cube adding a margine to the Beam Resolution cell to allow for 
% the rain drops to fall trough the through beam over time. The stretched 
% pulse will make the resolution cell larger and divided into range bins 
% according to specified input data. To get a resonable average the middle 
% range bin will serve as the "average" range bin 
% 

% Part II 

% From Rain Parameters: 

% Generating the fundamental T Matrix concisting of 
% - Diameter Vector and Number of drops in each diameter interval 
% using the Marshall-Palmer exponential approximation. 

% - the backscatter RCSs for each Diameter size 
% - the total number of drops for each size per m^3 

% Part III 

% Randomly place all drops in cube 
% Part IV 

% For every pulse the coherent Electric field return is calculated. Between 
% pulses all drops are moved. 

% Part V 

% Calculate and plot results 

% - 

clear all 
close all 
clc 

% Input Data 

c=3e8; % Speed of light [m/s] 

% Radar Data 

P_t= 10e3; % Transmitted power 

G_dB=20; % Maximum gain in dB 

G=10^(G_dB/10); % Maximum gain 

angle_3db=0.20; % 3 dB Beam angle [degrees], (half half power beamwidth) 

pulsewidth=0.20; % Pulsewidth [microsecond] 

elevation=3; % Elevation angle [degrees] 

margine=5; % margin for resolution cell[m] 

azimuth=0; % Azimuth angle [rad] 

PRT_l=l/200; % Pulse Repetition Time [seconds] 

distance=0.5; % Range from the radar to the resolution cell [Km] 
wavelength=0.1; % Wavelength [m] 

frequency=c/wavelength; % Frequency [Hz] 
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% - 

%Introduced in Pulse Compression version 
ehirp_width=2e6; % Frequeney step iir [Flz] 

ir_steps=10; % Number of Frequeney steps 

fq_step=ehirp_width/n_steps; % Frequeney step in [Flz] 

d_range=e*pulsewidth*n_steps*le-6/2/n_steps; % range bin due to pulse eompression 
num_mns=10; % Number of runs (pulse eompressed) 

range_bin=(distanee* 1000): drange: (distanee * 1000+(n_steps-1) *d_range); 

% - 

% Rain parameters 

diam_low=0.05; % Lowest diameter of Preeipitation [mm] 
diam_lim=l; % Divider between low and mid [mm] 

diam_mid=5; % Divider between mid and high [mm] 

diam_high=15; % Flighest diameter of Preeipitation [mm] 
resolution 1=0.08; % Lower diameter resolution [mm] 

resolution2=0.3; % Fligher diameter resolution [mm] 
re solutions =0.8; % Flighest diameter resolution [mm] 

sort_approx=l; % RCS approximation model; [l]=Rayleigh 
kw=0.93; % abs(K)*2 approx for water 

rain_rate=50; % Rain Rate [mm/h] 

rain_anglel=0; % degrees Rain Angle (Phi) 

rain_angle2=0; % degrees Rain Angle (Theta) 
wind_veetor=[-10;0;0]; % Wind speed veetor (x,y,z) [m/s] 
spread=l; % Wind speed spread [m/s], (all direetions) 

spread_2=4; % Drop veloeity spread [m/s] 

Norm=200; % Normalization eoeffieient 

% Direetional eosines 

u=sin(pi/2-elevation*pi/180)*eos(azimuth*pi/180); 
v=sin(pi/2-elevation*pi/180) * sin(azimuth*pi/l 80); 
w=eos(pi/2-elevation*pi/180); 

% - 

% Part I, Generate the resolution volume 

% - 

% Approximation of the resolution eell of the pulse eireular box, 

% for any azimuth angle 

% - 

% INPUT 
% - distanee 
% - angle SdB 
% - pulsewidth 
% - elevation 
% - margine 
% - azimuth 

% - 

% OUTPUT 

% - volume resolution: eomputed volume of the resolution eell 
% - volume box: volume of the box where the drops are put 
% - box: the 8 reetangular eoordinates of the box 
% - eoord vol res: the eoordinates of the resolution eell where the drops 
% are taken. 

% - 
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% - 

% Creating the eireular beam resolution eell 
range_res=(e*pulsewidth*n_steps*le-6)/2; 


% small angle approximation allows for following 
angle_rad=(angle_3 db *pi)/180; 

% 3 dB beamwidth [rad] 
radius 1=(tan(angle_rad)* (distanee *1000)); 

% Radius of beam at R=distanee 
radius2=(tan(angle_rad)*((distanee* 1000)+range_res)); ... 

% Radius of beam at R=distanee + range resolution 

surfaee 1 =pi *radius 1 "^2; 
surfaee2=pi *radius2*2; 

% volume of eone = base surfaee times height for all three 
volumeeone 1 ^(surfaee 1 *(distanee* 1000))/3; 
volume_eone2=(surfaee2*((distanee*1000)+range_res))/3; 
volume_resolution=(volume_eone2-volume_eonel); % Beam Cone Volume 

% To assure right volume for averaging 
DR=1:10; 

range_res_ref=(e*pulsewidth*DR* 1 e-6)/2; 

radius2_ref=(tan(angle_rad)*((distanee*1000)+range_res_ref)); ... 
surfaee2_ref=pi*radius2_ref^2; 

volume_eone2_ref=(surfaee2_ref.*((distanee*1000)+range_res_ref))/3; 

volume_resolution_ref=(volume_eone2_ref-[volume_eonel volume_eone2_ref(l:end-l)]); 

% Beam Cone Volume 


% eomputation of the resolution eell eoordinates, spherieal eoordinates 
% elevation in radians 

elevation_rad=(elevation*pi)/180; % Elevation [rad] 

azimuth_rad=(azimuth*pi)/180; % Azimuth [rad] 

% Matrix used when eomparing wheather a drop is inside or outside 
% beamwidth 

eoord_vol_res=[distanee*1000,(distanee*1000)+range_res;elevation_rad-angle_rad,... 

elevation_rad+angle_rad;azimuth_rad-angle_rad,azimuth_rad+angle_rad]; 

% - 


% - 

% Creating the box 

% eomputation of the box and its eoordinates 

eoinl=[azimuth_rad+angle_rad,elevation_rad+angle_rad,distanee* 1000]; 
eoin2=[azimuth_rad-angle_rad,elevation_rad+angle_rad,distanee* 1000]; 
eoin3=[azimuth_rad+angle_rad,elevation_rad-angle_rad,distanee *1000]; 
eoin4=[azimuth_rad-angle_rad,elevation_rad-angle_rad,distanee*1000]; 
eoin5=[azimuth_rad+angle_rad,elevation_rad+angle_rad,distanee* 1000+range_res]; 
eoin6=[azimuth_rad-angle_rad,elevation_rad+angle_rad,distanee*1000+range_res]; 
eoin7=[azimuth_rad+angle_rad,elevation_rad-angle_rad,distanee*1000+range_res]; 
eoin8=[azimuth_rad-angle_rad,elevation_rad-angle_rad,distanee* 1000+range_res]; 

[xl,yl,zl]=sph2eart(eoinl(l),eoinl(2),eoinl(3)); 

[x2,y2,z2]=sph2eart(eoin2(l),eoin2(2),eoin2(3)); 

141 






[x3 ,y3 ,z3]=sph2cart(coin3 (l),coin3 (2),coin3 (3)); 
[x4,y4,z4]=sph2cart(coin4(l),coin4(2),coin4(3)); 
[x5,y5,z5]=sph2cart(coin5(l),coin5(2),coin5(3)); 
[x6,y6,z6]=sph2cart(coin6(l),coin6(2),coin6(3)); 

[x7,y7,z7]=sph2cart(coin7 (l),coin7 (2),coin7 (3)); 
[x8,y8,z8]=sph2cart(coin8(l),coin8(2),coin8(3)); 

xmax=max([xl,x2,x3,x4,x5,x6,x7,x8]); 

xmin=miii([xl,x2,x3,x4,x5,x6,x7,x8]); 

ymax=max([yl,y2,y3,y4,y5,y6,y7,y8]); 

ymin=min([yl,y2,y3,y4,y5,y6,y7,y8]); 

zmax=max([zl,z2,z3,z4,z5,z6,z7,z8]); 

zmiii=miii([zl,z2,z3,z4,z5,z6,z7,z8]); 

% the box with its margin 

coin_l=[xmin-margine,ymin-margine,zmin-margine]; 

coin_2=[xmin-margine,ymax+margine,zmin-margine]; 

coin_3=[xmax+margine,ymax+margine,zmin-margine]; 

coin_4=[xmax+margine,ymin-margine,zmin-margine]; 

coin_5=[xmin-margine,ymin-margine,zmax+margine]; 

coin_6=[xmin-margine,ymax+margine,zmax+margine]; 

coin_7=[xmax+margine,ymax+margine,zmax+margine]; 

coin_8=[xmax+margine,ymin-margine,zmax+margine]; 

box=[coin r,coin_2',coin_3',coin_4',coin_5',coin_6',coin_7',coin_8']; 

% calculate the volume of the box 

volume_box=abs((xmax+margine-(xmin-margine))*(ymax+margine-(ymin-margine))*... 

(zmax+margine-(zmin-margine))); 


cube=[coinr,coin_2',coin_3',coin_4',coinr,coin_5',coin_6',coin_2',... 

coin_6',coin_7',coin_3',coin_7',coin_8',coin_4',coin_8',coin_5']; 

% - 

% - 

% Part II, Generate the fundamental T Matrix 

% - 

% This part produces a matrix of all drops that contribute to the total 
% reflectivity of the resolution cell 

% - 

% First all sizes are created, with two vectors of different size resolution. 

% The size vectors represents all diameters. Next the backscatter RCS for each 
% size is computed. 

% Using the Marshall-Palmer exponential approximation, the number of drops of 
% each diameter is computed. 

% - 

% INPUT 
% - diam low 
% - diam lim 
% - diam high 
% - resolution! 

% - resolution2 
% - sort approx 
% - wavelength 
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% - Kw 
% - rain_rate 

% - 

% OUTPUT 
% - number drops 
% Matrix T: 

% - first row all sizes 
% - second row all the backscatter RCSs 
% - third row the total number of drops for each size per m^3 
% the columns start with the smallest size 

% - 

% Create dropsize vector 

% Precipitaion Diameter vector (lowest<D<divider) 
vector_size_ 1 =diam_low resolution 1 :diam_lim; 

% Precipitaion Diameter vector (divider<D<mid) 
vector_size_2=diam_lim:resolution2:diam_mid; 

% Precipitaion Diameter vector (mid<D<highest) 
vector_size_3=diam_mid:resolution3:diam_high; 

% Precipitaion Diameter Vector 

vector_size=[vector_size_l,vector_size_2,vector_size_3]; 


% - 

% Create RCS matrix for all Rain drops of Precipitation Vector and for all 
% frequencies. 

fc|_rr=ff equency: fq_step: frequency+(n_steps-1 )* fq_step; 
rcs(l:n_steps,:)=((kw^2)*(pU5)/(cM)*fq_rr'.^4)*((vector_size*le-3).^6); 
%- 


% - 

% The Marshall Palmer approximation is used to create Precipitaion Diameter 
% Distribution Vector (Drop Size Distribution N(D) Vector) 

% The number of drops are scaled down by a factor (Norm) to improve the 
% speed of the simulator 

% - 

lambda=4.1 * (rain_rate^(-0.21)); % [mm-1 ] 

number=8000*exp(-lambda*vector_size)/Norm; 

vectorres 1 =ones( 1 ,length(vector_size_l ))*resolutionl; 
vector_res2=ones(l,length(vector_size_2))*resolution2; 
vector_res3 =ones( 1, length(vector_size_3)) *resolution3; 

% delta(Diameter) dD 

vector_res=[vector_resl,vector_res2,vector_res3]; 

% the approximate "true" number of drops N(D)dD vector 
number=number.*vector_res; 

% - 


%■ 


% create matrix T without res 
T=[ vectorsize; numb er]; 


%■ 
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% - 

% Change matrix 'T' to matrix 'matrix T Round' 

% - 

% This part computes the total number of drops by utilizing the resolution 
% cell and the matrix T Round 

% The only difference between the matrix T Round and T, is the fact that 
% matrix T Round has the total number of drops on the third row and they 
% are rounded of to nearest integer value. 

% - 


matrix_T_Round=T; 

matrix_T_Round(2, :)=round(matrix_T_Round(2, :)*volume_box); 

% Keep track of errors due to round off 
diff=T(2,:)*volume_box-matrix_T_Round(2,:); 

error_rcs=diff*rcs'*Norm/volume_box*sum(volume_resolution_ref)/n_steps; 

error_Z=error_rcs/sum(volume_resolution_ref)/n_steps*c''4./fc|_rr.^4/pi^5/... 

abs(kw)*2/le-18; 


% - 

% Part III, Initial positioning of the drops within the volume 

% - 

% This part generates the initial position of the drops in the volume of 
% the box. 

% The positions are recorded in the fde "sizexx.maf taking into account 
% the azimuth angle 

% - 

% INPUT 
% - box 

% - matrix T Round 
% - elevation 
% - azimuth 

% - 

% OUTPUT 

% All fdes including the sizes. Each contains the positions of all drops 
% for a given diameter. The data recoreded as follows: 

% - radial distance r 
% - elevation phi 
% - azimuth theta 

% - angle with respect to beam center. This parameter will be 
% used to identify drops which are in the beam. 

% - 


% Randomly put the calculated number of drops of each diameter in the box 
number_size=length(matrix_T_Round( 1,:; 

for i=l :number_size 

% select number of drops of certain Diameter 
drops=matrix_T_Round(2,i); 

% create as many positions on the x axis as # drops per size 
elements_x=rand( 1,drops); 

% put them in the range of interest 
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% Random placement x for every drop 

position_x=box( l,l)+(box(l,3)-box(l,l))* elementsx; 

% create as many positionson the y axis as # drops per size 
elements_y=rand( 1,drops); 

% put them in the range of interest 
% Random placement y for every drop 

position_y=box(2,l)+(box(2,2)-box(2,l))*elements_y; 

% create as many positions on the z axis as # drops per size 
elements_z=rand( 1,drops); 

% put them in the range of interest 
% Random placement z for every drop 

position_z=box(3,1 )+(box(3,5)-box(3,1 ))*elements_z; 

% matrix of results 

% Creating a drop Matrix using spherical coordinates 

[theta,phi,r]=cart2sph(position_x,position_y,position_z); 

% Calculating the range, phi, theta and angle to every drop with 
% reference to phi O and theta O. Used to estimate weighted power return. 
vector_azimuth=ones( 1 ,drops) * ((azimuth *pi)/180); 
vector_elevation=ones(l,drops)*((elevation*pi)/180); 
[xf,yf,zf|=sph2cart(vector_azimuth,vector_elevation,r); 
range=sqrt((position_x-xf). ^2 .+(position_y-y f) .^2. +(position_z-zf) .^2); 

% approximation 
angle=atan(range. /r); 

% Saving the drop positions in size%d.mat 
% matrixpos 

matrixpos=[r;phi;theta;angle]; 
save(sprintf('size%d.mat',i),'matrixpos'); 
end 

% - 


% - 

% Part IV, Calculate Coherent Electric field return. 
%- 


% For every pulse the coherent Electric field return is calculatet. Between 
% pulses all drops are moved. 


% 


% Run the simulation for specified number of pulses 
%- 


% Zero all vectors and matrices that are going to be used 
phase_vector=[]; 

sum_E=zeros(num_mns,n_steps); 

weighted_E=zeros(num_mns,n_steps); 

mnning_Z=[]; 

% Initiate reference variables 
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% Used to separate runs 


p=0; 

for runs=l :num_runs 
p=p+l; 

hh=0; % Used to separate frequeneies 

for fc|_new=frequeney: fq_step: frequeney+(n_steps-1) * fq_step 

hh=hh+l; 

k=2*pi*fq_new/e; 

disp(spriiitf('ealeulate the power-return of pulse %d\n',p)); 


% - 

% eomputes the (weighted) I and Q return 

% - 

% eomputes the total number of drops 

% - 

% INPUT 
% - Matrix T 
% - wavelength 
% - angle_3dB 
% - distanee 
% - eoord vol res 
% - number_pulses 

% - 

% OUTPUT 

% - total E: the reeeived eomplex Eleetrie field) 

% - power: the power 
% - voltage: the voltage 

% - drops beam: the number of drops in the beam 

% - 

% Only the drops whieh are in the resolution eell are seleeted for the 
% eomputation 

% - 


number_size=length(matrix_T_Round(l,:)); 

total_E=0; 

total_phase=0; 

new_E=0; 

for i=l :number_size 

% opening the file of registered data matrixpos. 
load(sprintf('size%d.mat',i)); 

% removing drops that are outside the resolution eell 
number_drops_per_size=length(matrixpos( 1,:)); 
dropmatrix=[]; 

% Filling the drop matrix with all drops inside the beam 

% resolution volume 

for d=l :number_drops_per_size 

if ((matrixpos(4,d)<=(angle_3db*pi/180))&(matrixpos(l,d)>=... 
range_bin(hh))&(matrixpos(l,d)<=range_bin(hh)-l-d_range)); 
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dropmatrix=[dropmatrix,matrixpos(:,d)]; 

end 

end 


% - 

% Convert to eartesian eoordinates 
if isempty(dropmatrix)==0 

[position_x,position_y,position_z]=sph2eart(dropmatrix(3, 
dropmatrix(2, :),dropmatrix( 1, 

% Create eomparison veetors 

veetor_azimuth=ones( 1, length(dropmatrix( 1,:))) * ((azimuth*pi)/180); 
veetor_elevation=ones( 1 ,length(dropmatrix( 1,: )))* ((elevation*pi)/180); 


% Create a phase veetor for the drops of eurrent Diameter 
r_phase=position_x. *u+position_y. *v+position_z. *w; 
phase=exp((2 *j *k). *r_phase); 

% Sum the eontributions to the E-field 

total_phase=total_phase+sum(phase); 

total_E=total_E+sum(sqrt(res(hh,i)).*phase); 


% - 

% eeate weighted power return, eoneidering range (r), theta, phi 
% with referenee to boresight and radar parameters 

E_weighted=sqrt(P_t)*G*e/fq_rr(hh)/(sqrt(4*pi)^3)./dropmatrix(l,:).^2.*... 
(exp(-4 *log(2). * ((dropmatrix(3,: )-veetor_azimuth) .^2/... 

(2 *angle_3db *pi/180)^2+((dropmatrix(2,: )- 
veetor_elevation) .^2/(2 *angle_3 db *pi/180)^2)))); 

new_E=new_E+sum(E_weighted. * (sqrt(res(hh,i)). *phase)); 

end 

end 


% - 

% add to E-field and phase veetor 
sum_E(p,hh)=total_E; 
weighted_E(p,hh)=new_E; 
phase_veetor(2*p-1 )=ANGLE(total_phase); 
%- 


end 

% - 

% This part simulates the movement of the drops and gives new positions. 
% First using the Doppler PRT (PRT l) 

% - 

% INPUT 

% - matrix T Round 
% - PRT 
% - azimuth 
% - elevation 
% - rain_anglel 
% - rain_angle2 
% - diam lim 
% - resolution! 
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% - resolution2 

% - wind vector 

% - 

% OUTPUTS 

% Output files containing sizes with new positions 
% structured the same way as in previous 

% there is also an output file for speed (one for each size) for the pulse 
% pair. 

% - 


numb er_size=length(matrix_T_Round( 1,:; 

% adding a Gaussian random term to make the wind speed spread 
if wind_vector( 1 )~=0 

new_wind_vector( 1 )=wind_vector( 1 )*normmd( 1 ,abs(spread/max(wind_vector( 1)))); 
else 

new_wind_vector( 1 )=0; 
end 

if wind_vector(2)~=0 

new_wind_vector(2)=wind_vector(2)*normmd(l,abs(spread/max(wind_vector(2)))); 

else 

new_wind_vector(2)=0; 
end 

if wind_vector(3)~=0 

new_wind_vector(3)=wind_vector(3)*normmd(l,abs(spread/max(wind_vector(3)))); 

else 

new_wind_vector(3)=0; 

end 

for i=l:number_size 

% opening of data file registered in matrixpos. 
load(sprintf('size%d.mat',i)); 

% total number of drops 
drops=matrix_T_Round(2,i); 

% transforming to rectangular coordinates 
[x,y,z]=sph2cart(matrixpos(3,:),matrixpos(2,:),matrixpos(l,:)); 

diameter=matrix_T_Round( 1 ,i); 


% - 

% Drop fall speed depends on drop diameter in the z axis, 

% Atlas-Ulbrich approximation is used 

% Approximation is valid in the diameter range 5*e-4, 5*e-3 
%- 


wtmax=386.6*((matrix_T_Round(l,i)*0.001)^0.67); 

velocityz=(wtmax*(ones(l,drops)))*cos((rain_anglel*pi)/180); 
velocityx=(wtmax* (ones( 1,drops))) * sin((rain_anglel*pi)/180)*... 
cos((rain_angle2*pi)/180); 

velocityy=(wtmax* (ones( 1,drops))) * sin((rain_anglel*pi)/180)*... 
sin((rain_angle2*pi)/l 80); 
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% - 

% - 

% 

% % adding another Gaussian random term to make the wind speed 

% % spread between different drops sizes 

% ifnew_wind_veetor(l)~=0 

% newest_wind_veetor( 1 )=new_wind_veetor( 1 )*normmd( 1 ,abs(spread. 

% /max(wind_veetor(l)))); 

% else 

% newest_wind_veetor(l)=0; 

% end 

% if new_wind_veetor(2)~=0 

% newest_wind_veetor(2)=new_wind_veetor(2)*normmd(l,abs(spread. 

% /max(wind_veetor(2)))); 

% else 

% newest_wind_veetor(2)=0; 

% end 

% if new_wind_vector(3)~=0 

% newest_wind_veetor(3)=new_wind_veetor(3)*normmd(l,abs(spread. 

% /max(wind_veetor(3)))); 

% else 

% newest_wind_veetor(3)=0; 

% end 

% 

% % If spread within pulse 

% veloeityxt=newest_wind_veetor( 1 )-veloeityx; 

% veloeityyt=newest_wind_veetor(2)-veloeityy; 

% veloeityzt=newest_wind_veetor(3)-veloeityz; 

% 

% 

% - 

% - 

% If spread for all drops regardless of size 
newest_wind_veetor=[]; 


% adding Gaussian random term to make the wind speed 
% spread between drops 
if new_wind_veetor(l)~=0 

newest_wind_veetor( 1, :)=new_wind_veetor( 1) *... 

normmd(l,abs(spread_2/max(new_wind_veetor(l))),l,drops); 

else 

newest_wind_veetor( 1,: )=0; 
end 

if new_wind_veetor(2)~=0 

newest_wind_veetor(2,:)=new_wind_veetor(2)*... 

normmd(l ,abs(spread_2/max(new_new_wind_veetor(2))), 1,drops); 

else 

newest_wind_veetor(2,: )=0; 
end 

if new_wind_veetor(3)~=0 

newest_wind_veetor(3, :)=new_wind_veetor(3) *... 

normmd(l,abs(spread_2/max(new_wind_veetor(3))),l,drops); 

else 

newest_wind_veetor(3,: )=0; 
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end 


% If spread within pulse 
veloeityxt=newest_wind_veetor(l,:)-veloeityx; 
veloeityyt=newest_wind_veetor(2,:)-veloeityy; 
veloeityzt=newest_wind_veetor(3,:)-veloeityz; 


%■ 

%■ 


% speed registration 

m_veloeity=[veloeityxt;veloeityyt;veloeityzt]; 

save(sprintf('veloeity%d.mat',i),'m_veloeity'); 

deplaeementx=velocityxt*PRT_l; 
deplaeementy=velocityyt*PRT_l; 
deplaeementz=veloeityzt*PRT_l; 

position_x=x+deplaeementx; 

position_y=y+deplaeementy; 

position_z=z+deplaeementz; 

% transformation into spherieal eoordinates 

[theta,phi,r]=eart2sph(position_x,position_y,position_z); 

% eomputation of new angles 
veetor_azimuth=ones( 1 ,drops) *(azimuth*pi)/180; 
veetor_elevation=ones(l,drops)*(elevation*pi)/180; 
[xf,yf,zf|=sph2eart(veetor_azimuth,veetor_elevation,r); 

range=sqrt((position_x-xf).^2.+(position_y-yf).^2.+(position_z-zf).^2); 
% approximation 
angle=atan(range ./r); 

% matrixpos 

matrixpos=[r;phi;theta;angle]; 

% save to file, matrixpos 

s=sprintf('size%d. mat',i); 
save(s, 'matrixpos'); 


end 

end 


%■ 


% - 

% ealeulate RCS, Z and Power return for every pulse ineluding error. The 
% refleetivity is eaeulated for eaeh volume bin and then summed and 
% averaged. 

% - 

for 0=1 :num_mns 

running_Z(o,:)=(Norm*abs(sum_E(o,:)).^2)*eM./fc|_rr(o)^4/pi^5/... 

abs(kw)^2/volume_resolution_ref(o)/le-18+error_Z(o); 

end 
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ranning_Z=sum(runiiing_Z,2)'/ii_steps; 


% Transform into a vector 


% - 

% Adding the error to the weighted E estimations 

% - 

error_power=error_rcs*P_t*G''2*c''2./fq_rr.^2/pi^2*(2*angle_3db*pi/l 80)^2*... 
(c*pulsewidth* le-6)/1024/log(2)/pi^2/(distance* le3)^2; 

for 0=1 :num_runs 

power_retum(o,:)=Norm*abs(weighted_E(o,:)).^2+error_power(o); 

end 

power_retum=sum(power_retum,2)'/n_steps; % Transform into a vector 

% to use to estimate zeroth moment 
for 0=1 :n_steps 

power_retum_ref_matrix(: ,o)=N orm*abs(weighted_E(: ,o)). ^2+error_power(o); 
end 

power_retum_ref=sum(power_retum_ref_matrix,l)/num_runs; % Transform into a vector 


% - 

% calculate average Power and Z 

% - 

Z_avg=sum(mnning_Z)/p; 

Z_avg_plot=ones( 1 ,p) *Z_avg; 

power_retum_avg=sum(power_retum)/p; % [W] 

power_retum_avg_plot=ones(l,p)*power_retum_avg; 

number_vector=l :p; 


% - 

% Calculate references 
% Integral form of Z 

% - 

ii=0:100/1000:100; 

zz=ii.^6.*8000.*exp(-4.1.*rain_rate *(-.21).*ii); 
real_Z=ones( l,p)*10*log 10(trapz(ii,zz)); 


% - 

% Evaluate average Power return to make Z estimation including error 

% - 

% To make the estimate with the correct frequency, the power return Matrix 
% is used. For the plot the result must me summed and divided by the number 
% of samples used to average (which means over both frequency and pulse). 

Z_estimation=power_retum_ref_matrix*((fq_rr.^2).*((c./fq_rr).^4))'/P_t/G''2/c^2*... 
pi^2/(2*angle_3 db *pi/180)^2/(c *pulsewidth* 1 e- 
6) * 1024*log(2) * (distance * 1 e3 +range_res/2)^2/... 

((kw''2)*(pC5))/le-18; 

Z_estimation_plot=ones(l,p)*10*logl0(sum(Z_estimation)/n_steps/num_mns); 


%■ 
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% Part V, Plot of results 
%- 


figure(l); 

kl=l:p; 

number_vector=kl.*PRT_l * 1000; 

text_string=['Power Return for all ',num2str(p),' pulse trains.'... 

Frequeney stepped from ',num2str(fc|_rr(l)/le9),' GHz to ',... 
num2str(fq_rr(n_steps)/le9),' GHz with ',num2str(fq_step/le6)... 

MHz steps.']; 

% also eonverting to dBm 

plot(number_veetor,10*logl0(power_retum)+30,'-xb', number_veetor,10*logl0... 

(power_retum_avg_plot)+30,'—r'); 
title(text_string); 
xlabel('Sample time [ms]'); 
ylabel('Power [dBm]'); 

legend('Power return', 'Average Power return',0) 
grid on 

figure (2) 

plot(number_veetor, 10*log 10(abs(running_Z)),'-xb',number_veetor, 10*... 
loglO(abs(Z_avg_plot)),'—r',number_veetor,real_Z,'—k',... 
numb er_veetor,Z_e stimation_p lot,' —b'); 
title('Estimated Z per pulse'); 
xlabel('Sample time'); 
ylabel('Z [dBZ]'); 

legend('Simulation Z, (instantaneous)', 'Estimated Z, (average retum)','Z=\int D^6N(D)dD',... 

'Estimated Z, Zeroth moment',0) 
grid on 
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